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Introduction. Despite the great amount of work that has been done on the solution 
of the partial differential equations of mathematical physics by separation of variables, 
the subject remains in an unfinished state. In a comprehensive treatment of the subject, 
the first problem is to determine separability conditions for the Laplace and Helmholtz 
equations and to find the coordinate systems that satisfy these conditions. The second 
problem is to tabulate all the separation equations. The first of these questions has been 
studied in recent papers [1], while the second is considered in the present contribution. 

Attempts have been made in the past to subsume all the separation equations under 
a single general equation; but these formulations have not been satisfactory, since they 
include equations that are never encountered in physics and they exclude equations 
that are definitely needed. For instance, the hypergeometric equation includes both 
Bessel and Legendre equations but excludes most of the separation equations that 
occur in field theory. There seems to be a widespread idea that Bécher’s “generalized 
Lamé equation’ with five singularities includes all the needed equations as confluent 
cases. According to Ince [2], “This systematization was suggested by the discovery of 
Klein and Bécher that the chief linear differential equations which arise out of the 
problems of mathematical physics can be derived from a single equation with five 
distinct regular singular points in which the difference between the two exponents 


relative to each singular point is }”’. 


Similar statements are made by Forsyth [3] and by Whittaker and Watson [4]. 
Actually, only a few of the required equations are obtainable from the Bécher prototype 
with five singularities (h = 5); though by employing a family of prototypes with h = 1 
to h = 7, one can obtain the necessary equations. 

The present study extends the general work of Bécher [5] by actually finding the 
separation equations for the eleven simply-separable systems of Eisenhart [6], the 
eleven symmetric cyclide systems, and eighteen cylindrical systems. All separation 
equations for these coordinate systems reduce to nineteen distinct equations of the Bécher 
type. 
Coordinate systems. Important partial differential equations of classical field 
theory are the Laplace equation, the Poisson equation, the diffusion equation, and the 
scalar wave equation. The solution of all these equations reduces essentially to the 
solution of a single partial differential equation, the Helmholtz equation: 

Viet «oe = 0. (1) 


*Received November 10, 1954. 
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For any orthogonal coordinate system (u’, u’, u*) in Euclidean 3-space, 














,_ifs 2, (2 2.) 2 
77" ” > du’ \g ;; du'] (2) 
If the Helmholtz equation is simply separable, 
g = Uv’) - U’) - UW), (3) 
and the three separation equations are 
1 ( a) > ” 
f;, du’ fi du* + dy ai; -s, (4) 


where a; are the separation constants and ®;,; are the elements of the Stackel matrix [1]. 

Simple separability of the Helmholtz equation occurs in the confocal quadric systems 

of Eisenhart [6]. Apparently, these eleven coordinate systems exhaust the possibilities 

for the Helmholtz equation. If x” = 0, however, we have the Laplace equation, which 

is R-separable [1] in a number of additional coordinates. The separation equations are 
again given by Eq. (4), but in place of Eq. (3) we have 
iF ee 72/,.2 73/..3 

aL: [a ee : (3a) 

Riu ,u,u) 








g=— 

The only R-separable coordinates that give promise of practical applications are 
those having considerable symmetry. They include toroidal coordinates, bispherical 
coordinates, and the four systems based on elliptic functions [7]. We include the rotation 
cyclides of Bécher [5] but deliberately omit his asymmetric cyclides as being too complex 
for use in physics and engineering. 

In this paper, therefore, we shall consider the separation equations for the following 
coordinate systems: 


I. Systems [8] allowing simple separation of the Helmholtz and Laplace equations. 


Cylindrical systems. 
1. Rectangular coordinates 
2. Circular-cylinder coordinates 
3. Elliptic-cylinder coordinates 
4. Parabolic-cylinder coordinates 


Rotational systems. 
5. Spherical coordinates 
6. Prolate spheroidal coordinates 
7. Oblate spheroidal coordinates 
8. Parabolic coordinates 


Asymmetric systems. 


9. Conical coordinates 
10. Ellipsoidal coordinates 
11. Paraboloidal coordinates 


II. Systems allowing R-separation of the Laplace equation. 


12. Tangent-sphere coordinates 
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13. Cardioid coordinates 

14. Bispherical coordinates 

15. Toroidal coordinates 

16. Inverse prolate spheroidal coordinates 
17. Inverse oblate spheroidal coordinates 
18. 6-Sphere coordinates 

19. Bi-cyclide coordinates (J1Rx)* 

20. Flat-ring cyclide coordinates (J1Ry) 
21. Disk-cyclide coordinates (J2R) 

22. Cap-cyclide coordinates (J3R) 


III. Cylindrical systems allowing simple separation of the Laplace equation in two 


dimensions only. 


4. Cardioid-cylinder coordinates (P3C) 

5. Hyperbolic-cylinder coordinates (P4C) 
26. Rose coordinates (P5C) 

27. Cassinian-oval coordinates (E2C) 

28. Inverse Cassinian-oval coordinates (E3C) 
29. Bi-cylindrical coordinates (E4C) 

30. Maxwell-cylinder coordinates (E5C) 

31. Logarithmic-cylinder coordinates (L1C) 
32. Ln tan cylinder coordinates (L2C) 

33. Ln cosh cylinder coordinates (L3C) 

34. Inverse elliptic-cylinder coordinates (H2C) 
35. sn-Cylinder coordinates (J1C) 

36. cn-Cylinder coordinates (J2C) 

37. Inverse sn-cylinder coordinates (J3C) 

38. Ln sn-cylinder coordinates (J4C) 

39. Ln en-cylinder coordinates (J5C) 

40. Zeta coordinates (J6C) 


Bécher equations. Except for one trivial case**, all the separation equations of 
field theory are second-order, linear differential equations that can be written in the 
form 

a’Z dZ ‘ - 
—3 + P(z) — + Q) Z = 0, (5) 
dz dz 


where P and Q are rational functions of z. We now introduce the restriction [10] that 





P(z) = Al ™ + = oe — Mn—1 |, 
2lLz—a, 2-@ Psgemae pe 
(6) 
I 
Q(z) =— Ay = Ay + + Az ; 
(Z-— ae 1s cS a,)”* bites (2 — i 


*See Reference [7 
tSee Reference [9]. 


**The equation in time obtained by separation of the diffusion equation. 
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where m; and / are positive integers. The equation has n singular points, which occur 
at 2 = 4; , 42, °*** G,-; , ©. The orders of the poles are m, , m, , --- m,_, ; and the 
order of the pole at © is equal to the greater of two integers giving the orders of 
(22 — 2°P(z)] and 2*Q(z) for z > ©. The total order of the poles ish = >°7., m,; . 

Equations defined by (5) and (6) are called Bécher equations. They are all obtainable 
from Bécher prototypes; but it seems simpler to define them directly by Eq. (6) than 
to employ the rather round-about method used by Bécher and Ince [2]. The direct 
approach also leads to a convenient method of classification [10] in terms of n instead 
of h. 

A Bécher equation may be specified by writing a sequence of integers representing 
the orders of the poles. The equation given by (5) and (6) is specified as 
{m, , M2, °** M,_, , m,}, where m, is the order of the pole at ~. The final integer m, 
holds a privileged position: it must not be moved; though the other integers may be 
interchanged at will. 

Some separation equations in their original form are not Bécher equations, but 
they can be made into Bécher equations by an elementary functional transformation 
of the independent variable. All equations obtained by such transformations may be 
regarded as equivalent, since if the solution of one is known, the solutions of the others 
are immediately evident. Among the equivalent equations obtained from a single equa- 
tion there is in general only one Bécher equation [10]. Thus the flexibility introduced 
by transformation of variables does not ordinarily cause any ambiguity in equation 
specification. 

The only exception to this rule occurs with very simple equations having 

(a) not more than two distinct singularities, 


(b) A, = A, = +: =A, =0. 


TABLE 1 


Classification of Bécher equations 
that are obtained by separation of variables in 40 coordinate systems 























n Name Original Degenerate Cases 
4 Heine equation* {1222} 
Wangerin equationt {1122} {1121} 
Lamé wave equation {1113} {1112}, {1111} 
3 Legendre wave equation {123} {122}, {121} 
Baer® wave equation {114} {113} 
2 Bessel wave equation {24} {23}, {22} 
Weber equation {14} 
Elementary {33} {32} 
1 Elementary {04} {01} 
Basic equations = 9, 


Total equations = 19. 





*E. Heine, Handbuch der Kugelfunctionen, Berlin, 1878, p. 445. 
tA. Wangerin, J. reine und angew. Math. vol. 82 (1877) p. 145. 
°K. Baer, Parabolische Koordinaten, Frankfurt, 1888. 
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TABLE 2 
The separation equations of field theory 


For each designation, the first line gives the Bécher equation, the second gives 
the transformation of variable that yields the equation of the third line. The 
fourth line presents the general solution. 


Differential equations and solutions 
dZ 


n Designation 
| 1 2 2 | 
+ + z— a;3_| dz 














aZ 
2+ z— a, 


{ 1222} 2 
arcane dz 





yp - & 





] 
2 
A 2 la 3 
+ | 4+ Aw + Az’ + Ay |z ot 


(2 — a,)(2 — a)"(z — as)” 
z = sn’ ft, 


UZ sng(dn®¢+ ken’ §) dZ 
ent dnt dt 


dy’ 
12 on? k’* sn* ¢ | = 
+ | 24 sn y 7 oe = cn t dn? t Z = 0, 





Heine functions. 

] l . dZ 
es 
z—-a@ 2-—4;j dz 


z—a, 


1 
2 
| Ao + Aye + Az2” kz = 0, 


(2 — a,)(2 — a2)(z — a)” 








_ _ VZ 
at: 


{1122} 
dz 





+ 





de? sn t ag 
+ E sn’ — a — aa( sn’ t+ = -) |z = (0, 


Wangerin functions. 

1 1 2 dZ 

5 91t 5 — 
bese det |. z2— a | 


| ; 
| Ain os Az |: asi 0, 


*lG-ele - ale - al 


Z2— a3 





aZ 





z = sn’ ft, 


~ 


CZ engdngdZ _ | 1 e 
ay?’ snt dt ut Z=0, 


Wangerin functions. 
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TABLE 2—(Continued) 








Differential equations and solutions 
































Designation 
a’Z 1 1 1 Z 
iia at 2 f+ +2 
Ao + Axe + Azz” _ 
* E — a,)(z — a,)(2 — 5 sali 
z= ¢", 
2 3853 __ a th ie 
(f° - OG oF +r 277 — (bh + oN] Ge 
+ [Pet —- ppt DP? + (0 +e)q]Z = 0, 
Lamé wave functions. 
a’Z 1 1 1 
imi det 2 i141 41-|¢ 
Ay + Az = 
” li — a,)(z2 — a,)(z — 5k =e 
a= “e 
2 2) / 4-2 2 a’Z a 2 dZ 
P= WY — oc) Gat Re — OF + ONG 
+ [0 +e)q — pp + DP)Z = 0, 
Z = AE;(¢) + BF;(S). 
| rie | 1 az 
‘ {1111} 4 i/_1_, 1,1 _|¢ 
Ao _ 
¥ li — a,)(z — a,)(z — ma —s 
= ¢, 
aZ adZ 





(° — BYE — eC) Fat $2"? — (8 +e) Te 


+ [(b? + c’)q]Z = 
Z = AEs) + BFS). 





It can be shown [10] that all these simple cases that occur in field theory reduce to the 
equations {04} and {01}. We stipulate that if a set of Bécher equations are equivalent, 
the specification of the set shall be that having smallest h and smallest n. 

Separation equations. A study of the separation equations for the 40 coordinate 
systems listed in this paper shows that there are only 19 distinct separation equations. 
The 19 include all the degenerate cases obtained by allowing one or more separation 
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TABLE 2—(Continued) 














n Designation Differential equations and solutions 
dz 1 } 1 2 | dZ 
3 {123} at 2 a—-  a= ae dz 





Ao + Axe + Az] ms 
+ [4 —ae-e)  ~® 


= 'Z 
1 @Z _ ab 
0-3-2 
+ [ wea — 2) +pp+1) - ole = 0, 
Legendre wave functions. 


aqZ.,1 1 2 dZ 
° an 2 +3l; + |2 


— a, Z2— a, 


+[-AtA oko. 


z2— a,)(z — a,)" 











2 
see; 


» @Z dZ 
(l— $) qa — 2S ae 





+[po+-72 |z=0, 











st 
Z = AP3(t) + BQ). 
a’Z 1 1 2 dZ 
. —_ #+3[-.+-2.]? 
Ae = 
+ Lait ak -o. 
z= °°, 
» ae dZ q° 
1- %)=5-2— - —+* 5-27 =0, 
(l— 9) G5 — 2 - Ge = 0 


Z = APs) + BQs(9). 





constants to be zero. They include also the equations obtained by separation of the 
vector Helmholtz equation [11] in the 5 coordinate systems in which this procedure 
has been shown to be possible. The distinct separation equations are designated in 
Table 1 and written in detail in Table 2. Note that by no means all of the equations 
have h = 5, as suggested by Ince. 

Some of the equations of Table 2 are familiar and their solutions are well-known. 
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TABLE 2—(Continued) 























n Designation Differential equations and solutions 
2 

3 {114} #2 41 L— 4 1_|z 
dz 2Lz2— a, zZz— a.| dz 


+ [Aet Ae + del " 


be 


$; 

ie — (b +0) dZ 

2(¢ b)(¢ — c) dt 

He - pp t+) - +097 


" (¢ — HS — oe) = 


Il 


Baer wave functions. 








az, i l 1 1Z 
. — di + 2 ; =a re =| z 


_Arst Az fo _ 
r E — a)(z — 5/2 ss 


z= cos ¢, 








i + (a — 2q¢ cos 26)Z = 0, 
dé 
Mathieu functions. 
2 9 7 j 2 
2 (24) Zi |_| dz | Ao t Ae + Awl wi 
dz 2 Lz — a,j dz (z — a) 
z= ’ 

2 r 2 

T4422 4 (i +¢- Piz = 0, 

di ¢ di 4 


Bessel wave functions. 





Others have been studied only slightly, if at all. In particular, separation of Laplace’s 
equation (x = 0) yields the Bessel equation, the Legendre equation, the Lamé equation. 
Separation of the Helmholtz equation (x # 0) yields analogous but somewhat more 
complicated equations which we have called the Bessel wave equation, the Legendre 
wave equation, the Lamé wave equation. Despite their importance, very little attention 
has been given to these equations. There are also two distinct equations obtained from 
rotational cyclide coordinates, and these equations likewise need further study. 
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TABLE 2—(Continued) 





Designation 





Differential equations and solutions 








fé i 
{23} 


bo 
bo 


, 14} 


130} 








i 2_|@ | As + Az), 2 
zy 1/2 a (z — a,)* 4=0, 








Z = AJ,(qz) + BY,(q). 











: dZ Me = 
2 | 2. | ds * | =: — [2 -— 


Z = Ag? + Br”. 


aZ 1 dZ | As + Ael, oa 


dz” 2(z — a,) dz z-a, 











Z ” . go 
ae t | a + 2) — “4-14 = 0, 


Z = AD,(qg) + BD_,-1(qS). 


= | 3_| aZ + | 42+ Ad ly = 0, 
z—a,Jj dz (z — a,) 





Z= 2 [A sin pt + B cos pf]. 


dz. if 3 ]azZ Ay 
de * 2 ; A dz * leer ” 











dZ,2dZ_ ppt+idy _ 
dy” + ¢ df e alien. 




















n Designation 
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TABLE 2—(Continued) 























Differential equations and solutions 
1 (04} we + AoZ = 0, 
z= ¢, 
aZ 26, 
de + pZ=0, 
Z = Asin pé + B cos pt. 
1 {01} FJ = 0, 
dz 
e= §; 
a’Z 
ag = 9 
Z=A+B 
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ON THE “BANG-BANG” CONTROL PROBLEM* 


BY 
R. BELLMAN, I. GLICKSBERG AND O. GROSS 
The RAND Corporation, Santa Monica, Calif. 


Summary. Let S be a physical system whose state at any time is described by an 
n-dimensional vector x(t), where x(t) is determined by a linear differential equation 
dz/dt = Az, with A a constant matrix. Application of external influences will yield an 
inhomogeneous equation, dz/dt = Az + f, where f, the “forcing term’’, represents the 
control. A problem of some importance in the theory of control circuits is that of choosing 
f so as to reduce z to 0 in minimum time. If f is restricted to belong to the class of vectors 
whose ith components can assume only the values + b,; , the control is said to be of the 
“bang-bang”’ type. 

Various aspects of the above problem have been treated by McDonald, Bushaw, 
LaSalle and Rose. We shall consider here the case where all the solutions of dz/dt = Az 
approach zero as t ©. In this case we prove that the problem of determining f so as 
to minimize the time required to transform the system into the rest position subject to 
the requirement that f, , the 7th component, satisfies the constraint | f; | < 6; may be 
reduced to the case where f,; = + b; . Furthermore, we show that if all the characteristic 
roots of A are real and negative, f; need change value only a finite number of times at 
most, dependent upon the dimension of the system. 

Finally, an example is given for n = 2, illustrating the procedure that can be followed 
and the results that can be obtained. 

1. Introduction. Let z be an n-dimensional vector function of ¢ satisfying the linear 
differential equation 
4 = Az+f, 2(0) =, (1.1) 
where we assume that: 

a. A isareal, constant matrix of order n, whose characteristic roots all have negative 

real parts; 

b. f is restricted to be real, measurable, and to have components satisfying the 

constraints, | f; | < 1. 

The first condition is the necessary and sufficient condition that all the solutions of 
(1.1) approach zero as t >. 

The problem we wish to consider is that of determining the vectors f which, subject 
to the constraint (b), reduce z to zero in minimum time. This is a problem of Bolza of 
rather unconventional type, and the techniques we shall employ are quite different 
from the classical ones. 

We shall establish two results: 


THEOREM 1. Under the above conditions, an f which reduces z to zero in minimum 
time exists, and has components f; for which | f; | = 1. 


*Received December 10, 1954; revised manuscript received March 14, 1955. 
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THEOREM 2. [f the characteristic roots of A are real, distinct, and negative, a minimizing 


f exists with components f; for which | f; | = 1, and each f; changes sign at most (n — 1) 
times. 


The statement in Theorem 1 has been assumed in the past on an intuitive basis, 
see McDonald, [3], and has been established in various cases by Bushaw [1], LaSalle 
[2], and Rose [4]. The only paper we have had access to is that by Rose, and his methods 
are distinct from ours. In addition, he is primarily interested in the case where the 
condition in (a) is not satisfied. 

Problems of this type arise in connection with many different types of control 
processes. A discussion of the connection with servomechanisms is sketched in [2]. 

2. Proof of Theorem 2. We shall consider in detail only the case of Theorem 2, 
where the characteristic roots of A are real and negative. It will be clear from the treat- 
ment of this case how the proof of Theorem 1 goes. 

Let X be a square matrix whose columns are the n linearly independent eigenvectors 
x; of A, and let \; (j = 1, --- , n) be the corresponding n distinct, negative eigenvalues 
of A; clearly, X is non-singular and all its elements are real. Finally, denote by A the 
diagonal matrix whose jth diagonal element is \; . We have 


Az; = i,2; , (2.1) 


whence we see that 4X = X A: hence 


If now in (2.1) we make the transformation z = Xy, we obtain using (2.2), 


y'(Q) = X c, 


(2.3) 
y(t) = Ay(t) + X'f (0), 
or, componentwise, 
y(t) = yA) + Do ai f(d, (2.4) 
7=1 
where the a’s are the elements of X~". Solving for y;(t), we obtain 
at n 
y(t) = y,(0) exp (A,;t) + exp (A;2) | exp (—\,s) S a;; f(s) ds. (2.5) 
J0 j=1 
Since z(t) = 0 is equivalent to y(t) = 0, we wish to find the least ¢ for which, for some 
f, yi = 0,7 = 1, --- , n,1.e., for which 
t n 
—y (0) = [ exp (—\,8) 2. a; ;f;(s) ds, t= 1,--- xn (2.6) 


for some f. 

Our first observation is that, given any starting value 4(0) ¥ 0 there exists at > 0 
and an f, such that (2.6) is satisfied. In fact, there is a constant vector f(s) = k which 
does the trick for some ¢ sufficiently large. For substituting f;(s) = k, in (2.6), we obtain 


y (0) a Ay(0) 


>” a; ;k; re > = 
2, mai — fo exp (—A,s) ds exp (—A,t) — 1’ 
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whence, by virtue of the definition of the a’s, 


Ay.(0) 
k: = Vii it ° 2.7 
F x exp (—A,é) — 1 —_ 
Since — A; > O the right member of (2.7) can be made as small in magnitude as we 
please for sufficiently large t, and hence we can insure that | k; | < 1. 
For each t > 0 we have a linear mapping p, taking f into the n-dimensional vector 
with 7th component 





[exp (-2) D eussil0) ds, (2.8) 


and this mapping clearly takes our basic convex set of f’s onto a convex subset C(t) of 
euclidean n-space. For any f in our basic set there is another, f, in the set which agrees 
with f for s < ¢ and vanishes for s > ¢, so that, for t’ > t, p.-f = pf = pf, by (2.8), 
and p,f is in C(t’). Thus C(@) increases with t. 

Now our desired least time is, by (2.6), the least ¢ > 0 for which C(é) contains the 
vector — »(0). Since C(t) increases, we have an interval (t) , ©) for which C(é) contains 
this vector, while for t < ¢, this is not the case. We can see that C(¢,) also contains this 
vector as follows. 

Denoting for any vector x = (a, , 2, «+: , %,) the euclidean norm (>>; 23)'” by 

x ||, we obtain, using (2.8), a constant k = k(t.) with the property that for all f and 
t, t' in a finite interval [0, t] we have || p,.f — pf || < k|t — v |; thus, for|¢t— ?’ | 
small every point of C(t’) is close to a point of C(t). Since — y(0) is in C(¢) forallt > %, 
— y(0) must be at zero distance from C(t.) so that if we show this set is closed — y(0) 
must actually be in it. But each C(d) is closed, since by a well known fact about Banach 
spaces [5], our basic set of f’s may be topologized so as to be compact and render each 
p, continuous. Thus C(é), as the continuous image of a compact set, is compact, hence 
closed. 

Let us return to the fact that — y(0) is not in C(t) for t < t . From the theory of 
convex sets [6] this implies that we have a vector 6‘ of unit norm, for which, in the usual 
inner product notation, (6', p,f) < [6@‘, — y(0)] for every f. Since the vectors of unit 
norm are compact in the euclidean topology, we may select a sequence ¢, increasing 
to t, for which 6" converges to some vector @ of unit norm. But since p,,f converges to 
p.,.f, (0, pf) lim (6°", pf) < lim [6°, — y(0)] = [6, — y(0)]. Thus if f* denotes an 


f for which p,,f* = — y(0) we have (6, p.,f) < (6, p:,f*) for all f, hence constants 6, , 
-- , @, , not all zero for which f* maximizes the expression 
» 4, | exp (—A,8) >> a; f(s) ds = >> [ (>> 0.a;; exp [—A,8]) f(s) ds. (2.9) 


But this expression clearly has as its maximum 


> [ |S @,a;; exp (—d,8) | ds (2.10) 


achieved by setting f;(s) = sgn } 6,a;, exp [— A,s]). Thus it is clear that f¥(s) = 
sgn CG. 6,a;; exp |— A,s]) almost everywhere on the set where » a 6,a,;; exp (— A,s) + 0. 

Our principal result now follows, namely that we can achieve minimal time by 
restricting f to assume componentwise + 1 on a finite number of intervals; in fact, in 
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the case considered, each component need change sign at most n — 1 times. This latter 
statement is a simple consequence of the fact that unless the continuous function ¢; 
given by ¢,;(s) = >-"-, 6,a;; exp (— ,8) is identically zero (in which case it makes no 
difference as to our choice of f*), it can have at most n — 1 real zeros. This is well known 
and there is a simple inductive proof. 

3. A special case of n = 2. Consider the problem as before, with 


—3 —2 
A= 
l 0 
thus, 
ei = —32z, —-22+f:, 2422+ f2. (3.1) 
The transformation 
2=YM—- ye, 2= —"W + Yy (3.2) 
reduces the above system to 
y= —Mthtih, w=—wetht2f, (3.3) 


and we obtain, as before, for the set of admissible starting values, for a given tand f,, fo, 
at 


—y,(0) = | ef) + f.9)] ds, 
(3.4) 


at 


—y,(0) = | e'[fuls) + 2f20)] ds. 


JO 

From the preceding section, we know that if ¢* is minimal, then the optimal f* is given by 
f.(s) = sgn (6,e°" + 6.€°), 

(3.5) 
f.(s) = sgn (6,e°°+ 26.¢’). 


If we now ask the question ‘For what set of starting values y is it optimal to choose 
f: = 1, fe = 1 on an initial interval?” with a similar question for the other combinations 
+ 1, it is readily seen that the answers will determine an optimal policy. This is clear, 
since any continuation of an optimal policy must be again optimal with respect to the 
new starting values. We thus have 
ase 
—y,(0) = | e** {sgn (0,e°" + 62e°) + sgn (6,e°" + 26,e")} ds, 


“0 


(3.6) 
—y(0) = / e*{sgn (0,e7" + 6,e°) + 2 sgn (6,e°° + 20,€")} ds. 


0 


To answer the first question, for what values of y is it optimal to set f, = f, = 1 on an 
initial interval, we note that this is equivalent to the conditions 


i* > 0, 
6, + 0, > 0, (3.7) 
6, + 20, > 0. 
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Now, since the functions 6,e"* + 6,e°, 0,e"* + 26,e° can each vanish at most once, 
we see that the above case breaks down into four sub-cases, namely: 
(a) 6, exp (2¢*) + 6, exp (t*) > 0, oe ws Se (ce) <,>, oS << 
(3.8) 
6, exp (2/*) + 26, exp (t*) > 0. 
Case (a) is trivial and consists of the arc a’ illustrated in Fig. 1. a’ is defined para- 
metrically by 
yi(0) = 1 — exp (2¢*) 
i*>0 (3.9) 
y2(0) = 3(1 — exp [é*)), 





¥,(0) 
43 








f=! 
fee! 

















Fie. 1. 


as one can readily verify by working out the integrals. Moreover, the curve defines an 
optimal path, since the solution of the differential equation is, with f,; = f. = 1 identi- 
cally, and y,(0), y2(0) defined as above, precisely a sub-are of a’ beginning at y(0) and 
terminating at the origin. 
Case (b) is vacuous, for if we have 
6, exp (2¢*) + 6, exp (t*) > 0 and 
(3.10) 

6, exp (2¢*) + 26, exp (i*) < 0, 
we obtain, by subtraction, and the condition ¢* > 0, that 6, < 0. But, 6, + @ > 0, 
whence 6, > 0. We thus have 


6, exp (2t*) + 26, exp (t*) > 6, exp (¢*) + 26, exp (¢*) 
= exp (/*) (6, + 26.) >0 (8.11) 
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Z,(0) 
43 








f=! 
1,21 














which contradicts 
6, exp (2t*) + 26, exp (t*) < 0. (3.12) 
We shall treat case (3.8c) in detail. Case (3.8d) can be treated similarly, but is a 
trifle more involved, albeit elementary, and will be omitted on those grounds. 


We have, upon substituting in (3.6) for case (c): 


ain . at® 
—y,(0) = | exp (2s) ds — | exp (2s) ds + | exp (2s) ds, 
70 Jin 64/64 0 
(3.13) 
In at nt® 
—y(0) = | exp (s) ds — | exp (s) ds + 2 | exp (s) ds. 
70 /in(-—86 ) 70 


Simplifying, we obtain 
—y,(0) = (6, 0,)° ~ i, 


(3.14) 
—y(0) = —2(6./6,) + exp (¢*) — 3. 
If now we set x* = exp (¢*), our conditions become 
ee 
6+ 6> 0, 
6, + 26, > 0, 
(3.15) 


6,2* + 0 < 0, 

6,2* + 20, > 0, 
yi(0) = 1 — (6,/6,)", 
y2(0) = 3 + 2(0./0,) — x*. 
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We easily obtain from the above that @, < 0. By homogeneity, we can set @, = — 1, 
6, = \ and we obtain the equivalent conditions 


A>sz*>a>I1 (A) 
yi(0) = 1 — 2’ 
(B) 
y(0) = 3 — 2A — a* 


i.e., we wish to find the image of all pairs (x*, A) satisfying (A) under the"mapping 
defined by (B). Pictorially this is represented by Fig. 3. 








N Y aaas — 
Lk. — 















































| 
| 
l 
{ | x* 
2 


On the other hand, the Jacobian of the transformation (B) is given by 


—2r 0 
J, = = 2. ~ 0 (3.16) 
«if | 


throughout (A); hence the transformation is non-singular and the boundary of the 
image is the image of the boundary. Making use of this fact we obtain the region for 
case (3.8¢): 

y:i(0) <0 (3.17) 


and 
3 — 4[1 — y,(0)]'” < y.0) < 3 — 3[1 — y,)]’”. (3.18) 


In a similar manner we obtain a region for case (3.8d). The union of cases (3.8a) through 
(3.8d) is the set of all starting values for which f, = f. = 1 is optimal on an initial 
interval. In a similar manner we obtain the region f; = 1, f. = — 1. (Notice that we 
need not compute the other regions since they can be obtained by skew-symmetry.) 

The final result of our calculations is illustrated in Figs. 1 and 2. Figure 2 is the 
image of Fig. 1 under our initial transformation and gives the optimal policy in terms 
of our initial starting vector c = [z,(0), z.(0)]. 

In terms of optimal paths (see Fig. 2) we can state the following: A path initiating 
in the (1,1) region continues with f; = 1, f. = 1 until it strikes either the straight 
segment OB or the parabolic are 8. In the former case f, switches from 1 to — 1 and the 
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path continues along OB to the origin. In the latter case f; switches to — 1 at 8 and the 
path continues in the (— 1, 1) region until it intercepts the parabolic arc a at which 
f, changes from 1 to —1 and a is followed to the origin with f, = f, = — 1. Similar 
remarks hold for the skew-symmetric regions. 
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THE FREE OSCILLATIONS OF A BUCKLED PANEL* 


BY 
W. D. HAYES 
Princeton University 
AND 
J. W. MILES 


University of California** 


Summary. The non-linear equations of motion of a buckled, two-dimensional 
panel are formulated in dimensionless form. A Fourier expansion in the space variable 
is introduced and an approximate solution is obtained for the period of free oscillation 
on the assumption of only two degrees of freedom, one of which is eliminated by the 
buckling constraint. The periods of the simplest asymmetric and symmetric modes are 
plotted as a function of the energy level. It is found that for very small energy levels 
the period of the buckled panel lies between the periods of the (two) unbuckled degrees 
of freedom. As the energy is increased the period approaches infinity at some critical 
energy level and thereafter decreases monotonically. 

1. Introduction. The subject of this study is the oscillation of a two-dimensional, 
buckled plate. The problem is of considerable practical interest and serves as a simple 
model of the phenomenon often described as “oil canning.’ Moreover, it is of some theo- 
retical interest that the original partial differential equation of motion is linear and that 
the deflection is required to be small but that the problem is nevertheless non-linear in 
consequence of the buckling constraint. 

The rather more difficult problem of the aerodynamic excitation of a buckled plate 
has been studied previously’, but the simplicity of the present problem allows a more 
complete solution. In particular, only the free oscillations will be treated, whereas 
self-excited oscillations are possible under aerodynamic excitation. 

2. Equations of motion. Consider a two-dimensional plate of flexural stiffness 
D and mass m per unit area that is loaded axially by two equal and opposite compressive 
forces F(t) per unit width, applied at the pin supports x = 0 and x = I. The coordinate 
x is measured along the plate from the left support, while y, the displacement, is measured 
positive up, as shown in Fig. 1; ¢ is the time. The equation of motion of the plate then 


reads 
a” 2s) 7 Hy ay _ 
Ox” (v Ox® +? Ox” sie _ 6. (2.1) 
y 
F( vais: aces 
me * L i 





Fig. 1. Two-dimensional buckled plate. 
*Received January 31, 1955. 
**The present work was carried out while the latter author was a Consultant to the Douglas Air- 
craft Co., Santa Monica, Calif., and appeared as Douglas Report SM-18439. 
iW. D. Hayes, A buckled plate in a supersonic stream, No. Amer. Avia. Rep. AL-1029, Downey, 


Calif., 1950. 
2J. W. Miles, Dynamic chordwise stability, No. Amer. Avia. Rep. AL-1140, Downey, Calif., 1950. 
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In general, D and m may be functions of x, but in the following they will be assumed 


to be constant. 
It is required to obtain a solution to (2.1) subject to the boundary conditions of 


zero displacement and zero bending moment at the ends, viz. 


y(0) = y(2) = 0, (2.2a) 
y.2(0) = y,,(l) = 0, (2.2b) 


and, in addition, to the constraint that the plate be buckled such that its length has the 
constant value / + él. It will be assumed that 4, the fraction of the length by which the 
plate is buckled, is small, whence the deflections also must be small; then, by definition, 


seat ] al au)’ 1/2 4 1 l (2x) oe 
head | E a (2u dx —l|= yi [ a dx. (2.3) 


In order to render the formulation dimensionless, the following dimensionless quan- 
tities and characteristic measures are introduced: 


= 2/l,, (2.4) 
n = y/lo = n€, 7), (2.5) 
t= t/t, (2.6) 
f() = F()/Fo , (2.7) 
lo = I/z, (2.8) 
to = (1/n)*(m/D)"”, (2.9) 
Fy = (n/l)’D. (2.10) 


The reference length, /, , is chosen in anticipation of the Fourier expansion in ¢é intro- 
duced below (Sec. 3). The reference time, ¢) , is the reciprocal of the fundamental (circular) 
frequency of simple harmonic motion for the plate in the absence of the buckling loads. 
F, is the Euler load, i.e., the smallest axial load required to produce static buckling in 
the absence of normal loading. 

Introducing the dimensionless variables ~, n, 7, and f in the original formulation 
of (1)-(3) leads to 


d*n an , On 
ge TIM eet ae = 9 1 = a, 9), (2.11) 
(0, r) = n(x, r) = O, (2.12a) 
nee(O, 7) = nee(r, 7) = O, (2.12b) 

1 f* oa’ 
= | (22 dg = 6. (2.13) 


3. Fourier series reduction. In order to eliminate the space variable from (2.11) 
and (2.13) it is convenient to expand 7(é) in a set of functions that individually satisfy 
the boundary conditions (2.12). (Such an expansion is permissible in virtue of the 
linearity of (2.11) in 7 as a function of £.) It is evident that the required set is furnished by 


naA(é) = sin (né), nm=mil,2,---. (3.1) 
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Accordingly, the Fourier expansion 


N 


n(t, 7) = > a,(7) sin (né) (3.2) 


n=1 
is posed. 
Introducing (2) in (2.11) and (2.13) and denoting differentiation with respect to r 
by dots yields 
a, (r) + [n* — n*f(r)Ja,(7) = 0, n=1,2,---,N (3.3) 
N 


> n’a 


1 


\— 


aw 


= 6 (3.4) 


_— 


which constitute a set of N + 1 equations for the determination of the N + 1 unknowns 
@, , Ge, °°** Qn, f. 
The initial values of » and n, are 


N 


n(é, 0) = >> a,(0) sin (né), (3.5a) 
n.(é, 0) = >< a;(0) sin (ng). (3.5b) 


In principle, N must be infinite for arbitrarily prescribed initial values of » and 7, , 
but, whatever the finite value of N that will permit a satisfactory approximation to 
these initial values, they must satisfy (3.4) and its time derivative, viz. 

N 


n’ar(0) = 46, (3.6a) 


1 


> n*a,(0)a,(0) = 0. (3.6b) 


It is clear that the minimum number of modes required to satisfy (6a) and (6b) in a 
non-trivial way is 2, a single mode leading only to a purely static solution. 

It does not appear possible to obtain an explicit, general solution to the problem 
formulated in (3)-(5). Accordingly, the analysis will be restricted to the simplest possible 
situation, where only two (Fourier) modes are required, say n = 1 and n = N(N > 2). 
For definiteness it will be assumed that ay(0) = 0 (although this limitation can be 
removed simply by changing the origin of time) and the initial conditions 


a,(0) = 28°”, ay(0) = 0, (3.7a) 
a,(0) = 0, ay(0) = 2(6w)'”’, (3.7b) 


posed, where, as will be brought out in more detail below, w defines the energy level of 
the oscillation. It should be emphasized that this reduced two degree of freedom problem 
does not afford a general solution inasmuch as it is not possible to superpose solutions 
for different NV. 

4. Energies. Before attempting an approximate solution to the problem formulated 
in the preceding section, it is instructive to consider the various energies involved. 
Inasmuch as the buckled state of the plate stands as a central feature of the analysis, 
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an appropriate reference quantity or unit of measurement is the work required to buckle 
the plate in the absence of normal loading, namely 


2 
W, = (=P). (4.1) 
The potential and kinetic energies per unit width of the plate are given by 

13 au) i 
K = 41 m( 2 dz, (4.2) 
’ 1 al “yy 

==> a (4 € 
V=5 | p(y dee. (4.3) 


Introducing the dimensionless notation of (2.4)-(2.10) and (4.1), these energies can be 


placed in the dimensionless forms 


ee ] * (an\° 
k = K/W. = = [ (2 dé, (4.4) 
276 Jo \Or 
P x 1 ar On 2 
v= V/Wo=5—| (<3) at. (4.5 
2rb Jo \ OE : ) 
Introducing the Fourier expansion of (3.2) in (4) and (5) yields 
1 N 
k(r) = — > [a;(7) I)’, (4.6) 
46 1 
i< 
v(r) = — Do n‘a2(z). (4.7) 
464 
Imposing the requirement of conservation of energy, viz. 
k+v=C (4.8) 
then yields 
vr. 1 4,2) Y f 
,D (a, +na,) = C. (4.9) 
1 


The constant C is determined by the initial energy in the system. For the two degree 
of freedom system having the initial conditions (3.7) 
(a; + ay) + (ai + Nas) = 48(1 + w), (4.10) 


where w is defined as the (dimensionless) energy in excess of the work required to produce 
the static configuration. 

Alternatively, (9) can be obtained by multiplying both sides of (3.3) by a, , summing 
over n, imposing the constraint that the derivative of (3.4) with respect to 7 vanish 
identically, and integrating the remainder with respect to r. 

5. Two mode solution. Returning to the solution of the two mode problem having 
the initial conditions of (3.7), the resulting equations of motion then can be obtained 
either from (3.3) and (3.4) or (4.10) and (3.4). The latter automatically eliminates the 
unknown f(r), which, however, can be determined from (3.3) once a,(r) and ay(7) are 


known. 
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The buckling constraint (3.4) will be satisfied by writing 
a,(r) = 26'” cos [A(7)], (5.1a) 
Nay(r) = 26'” sin [6(7)]. (5.1b) 


Il 


Substituting these expressions in (4.10) yields 
[1 + (N® — 1) sin’ 6][1 + (@°/N)*] — 1 = w. (5.2) 
It may be remarked that, by definition, w vanishes for the static configuration (@ = 0). 
In order to determine the period of the free oscillation of the plate as a function of 

its energy level, (2) is first solved for 0°, viz. 


We! | w — (N* — 1) sin’ 6 |” ‘ 
ian =n| * + (N? — 1)sin’ | ‘ (5.3) 





It is evident that the phase plane (@, 6) trajectories are closed, quasielliptic, curves 
for w < (N*® — 1) and single valued, periodic curves for w > (N* — 1). The critical 
trajectory w = (N° — 1) is designated as the “separatrix.” (The terminology is that of 
Minorsky”.) These topological distinctions are intimately connected with the singularities 
of the differential equation satisfied by the phase plane trajectories, namely [from 
logarithmic differentiation of (5.2)] 
6 dd ain aes a —sin 6 cos 6d 5.4 
(N?— DIN? + 0) (LF (= Din” ~— 





It is seen that the singular points (at which d6°/dé@ is indeterminate) arise at 6 = sx/2 
where s is any integer. Even values of s correspond to saddle points (points of stable, 
static equilibrium), in the neighborhood of which the trajectories defined by (3) are 
elliptic and w < 1. Odd values of s correspond to vortex points (points of unstable 
static equilibrium), in the neighborhood of which the trajectories are hyperbolic and 
| w— (N*® — 1) | «1. The separatrix passes through these vortex points and in this 
neighborhood coincides with the asymptotes of the hyperbolic trajectories. 


6 











Fic. 2. Phase plane trajectories of free oscillations (V = 2). 


8N. Minorsky, Non-linear mechanics, Edwards Bros., Ann Arbor, 1947. 
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Three typical trajectories obtained by assigning the values w = 3/2, 3 and 6 and 
N = 2 are plotted in Fig. 2. The pattern is, of course, repeated in each 7 interval of 6. 

The three possible types of oscillations under consideration are somewhat analogous 
to the “wig-wag”’ oscillations of a pendulum (total energy < 2mgl), the “loop the loop” 
oscillations (total energy > 2mgl) and the critical case where the total energy (= 2mgl) 
is just sufficient to invert the pendulum. By analogy with this classical problem, we 


may expect the period of oscillation of the buckled plate to increase from its initial 


value at w = 0 to infinity at w = N’ — 1 and then to decrease monotonically as w is 
increased beyond the critical value. As in the case of the pendulum, the maximum possible 
potential energy of the system (v = N’, of which 1 is already present in consequence 


of the buckling) is attained for w = N* — 1, and additional energy can appear only in 
kinetic form, whence the above inference of a monotonic decrease in the period. 

It remains to integrate (3). If w < N’ — 1 the period is defined as the time required 
to cover the complete orbit in the phase plane, which, in virtue of symmetry, is given by 


}i/s 


sin~*[{w/(N?-1)] 
T, =4 [ =. 


— w<N’—1. (5.5a) 
0 6°(6) 


If w > N’ — 1, the time required for a given configuration to be reestablished is evidently 
T=4f “, w>N*-1. (5.5b) 
In the first case it is convenient to introduce the change of variable 
1/2 
sin 9 = (= _ sin a (5.6) 
nt] i 


which, after substituting 6° from (3), reduces (5a) to 


r= 4] | prt tein a] a 5.7 
0 


N J (N? — 1) — w sin? a 


on 
“J 
6 
9 

Ww 


The direct substitution of 6° in (5b) yields 


. GPT. oe =. ee sy . aonite 
= — +, * (5.7 
Tr N Jo E — (N° — 1) sin’ 6 sd 5.7) 





Following standard procedures, each of the integrals in (7) reduces to an elliptic 
integral of the third kind, the end results being 


—— 4 ; tay" EF N a (.¥ ) 2 fs wy] _— 
‘iy = N (4, ve 1 I (N? a age 1 + Ww ’ N? = 1 ’ (5.8a) 
72 1/2 1/2 re 1/2 
au | (its) (A —1)"|, (5.8b) 
N Ww Ww 


. gp - [ Ga#sin’ py! 56 
I(k, B) iad i (1 a B° sin’ ¢) dy. (5.9) 


= 
L~ 
| 
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(It is of interest to remark that the arguments of J in (8a) and (8b) are the reciprocals 
of one another.) The integral J is evaluated as* 


.2 st 2\ 1/2 
I(k, 8) = K(k) + ; (¥=4) {E(k)F[k, sin’ (8/k)] — K(kK)E[k, sin“ (8/k)]}. (5.10) 
If w approximates the critical value (N* — 1) (5.10) is indeterminate, but a direct reduc- 
tion of either of (5.8) yields 


e —_— _4(N? — 1) _| ‘ , 
T = wt In Per os + O[w — (N° — 1]. (5.11) 
In the limiting cases of very small or very large values of w Eqs. (7) reduce to 
. = ae , 1 ) ) 2 Fr 19 
T. = nN? — Wg E +3 (<i —¢ | + Ow’), (5.12a) 
T, ite ag] "=D Nn + O(w*””). (5.12b) 


Now, in the absence of buckling, the periods of the first and Nth modes of natural oscilla- 
tion would be [set 6 = f = b, = 0 in (3.3) and (3.4)] 27 and 27/N’, and T,, lies between 
these values in the limit w = 0, although closer to the second than the first. This rather 
curious state of affairs, which, however, is typical of non-conservative systems, is to be 
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| Fic. 3. The period (7) of free oscillation of a buckled plate, referred to the fundamental period 
































(7) of the unbuckled plate. 


‘D. Bierens de Haan, Nouvelles tables d’integrales définies, G. E. Stechert & Co., New York, 1939, 


Sec. 57(1) and Sec. 61(8). 








26 W. D. HAYES AND J. W. MILES [Vol. XIV, No. 1 
contrasted with the coupled oscillations of a linear two degree of freedom system, where 
the resulting period cannot be between the periods of the uncoupled oscillations. 

The period obtained by dividing (5.8)-(5.12) by 2m (so that unity corresponds to the 
fundamental period of the unbuckled plate) is plotted vs. w in Fig. 3 for the two most 
important possibilities N = 2 (simplest asymmetric) and N = 3 (simplest symmetric). 
The reciprocal of this period, being the ratio of the frequency of oscillation to that for 
the unbuckled plate, is plotted in Fig. 4. 
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Fig. 4. The frequency (f) of free oscillation of a buckled plate, referred to'the fundamental frequency 
(f:) of the unbuckled plate. 












































THE STABILITY OF VISCOUS FLOW BETWEEN ROTATING CYLINDERS* 


BY 
H. STEINMAN 
Yerkes Observatory, University of Chicago 


Abstract. In this paper Chandrasekhar’s method of solving the characteristic 
value problem in the theory of the stability of viscous flow between rotating cylinders is 
carried to a higher approximation. A comparison between the results of the various 
calculations among themselves and with the experiments confirms the greater feasibility 
of the present method. 

I. Introduction. Recently S. Chandrasekhar (1954) [1] has investigated the problem 
of the stability of viscous flow between two concentric rotating cylinders. He employed a 
modified and more general method than had previously been used. The problem was 
first studied extensively by G. I. Taylor (1923). Later D. Meksyn (1946) expanded the 
theoretical treatment to cover a different range of values for the ratio of angular velocities. 
The significance of Chandrasekhar’s method lies in its relative simplicity and its ability 
to treat the widest range of values for the relevant parameters. This paper is an extension 
of Chandrasekhar’s work to obtain the solution of the problem correct to the second 
order in y/d, where y is the distance of a point from the inner cylinder and d is the differ- 
ence in the radii of the cylinders. 

The stationary solution of the hydrodynamical equations appropriate to the problem 


on hand is 
Vo) = ar +2, (1) 


where V(r) is the rotational velocity at a distance r from the axis of rotation; A and B 
are two constants related to the angular velocities 2, and 2, with which the cylinders 


of radii R, and R, are rotated by 


_ o L—#R2/Ri _ Ril —w) 
A= Q, 1 — R2/R? and B= Q, 1 — RR?’ (2) 








and 
si = 7a, (3) 


is the ratio of the angular velocities. 

By considering a symmetric perturbation of solution (1) by a periodic disturbance 
in the direction parallel to the axis of rotation with a wave length X, a characteristic 
value problem in a differential equation of order six is obtained (see Chandrasekhar, 
loc. cit., Eq. (4)). When the difference in radii of the two cylinders is small compared to 
their mean, we may, to the first order in y/d, replace A + B/r’ by 


A+ B/r’ = a1 — (1 — wel, (4) 
where 
r—R, 
t=y/d= R,— R, (5) 


*Received February 3, 1955. 
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This is the approximation on which both Meksyn and Chandrasekhar have solved the 
characteristic value problem. To the second order in y/d we have instead 


A+ B/r = o,f —(1l- uy) + : we Ro, + : (1 — p) a |: (6) 
by} F 1 


and we shall write this expression in the form 
A+ B/r’ = Q,[1 + at + BE’), (7) 


where 


3 (R. — R;) 
a= (u _~ pf + = oa Ba | 


3 Oo# R,) 
B = ~ p= (8) 
whe=1l+artéB. 


With A + B/r’ given by Eq. (6) the differential equation governing the problem 
is (see Chandrasekhar, Eq. (10)) 


(D? — a’)*v = —a’'T(1 + aé + BED, (9) 
where 
D = d/d, h = a(R, — R,) 
and 
O.d* 
p= —149¢ (10) 
y 


is the Taylor number. (In Eq. (10) » is the kinematic viscosity.) The six boundary 
conditions with respect to which Eq. (9) must be solved are 

v= (D’ — av = D(D’ -— av = 0 at gé= 0 and = 1. 11) 
Equation (9) together with the boundary conditions (11) determine 7 for any assigned a’. 
The critical Taylor number 7’. , at which instability sets in is then the minimum value 
of T as a function of a’. 

II. The revised method. Chandrasekhar’s method of solving the basic characteristic 
value problem consists in representing Eq. (9) as a pair of differential equations, one of 
order four and the other of order two and in arranging that the boundary conditions 
split in the same way; the fourth order equation together with the boundary conditions 
which go with it is then solved exactly. Thus, with the same substitutions as Chandra- 
sekhar (loc. cit., Eqs. (16)-(24)), we reduce the characteristic value problem to solving 
the pair of equations 

(D? — a’)’W = (1 + aé + BE) (12) 
and 
(D* — a)y¥ = —aTW, (13) 
together with the boundary conditions 


W DW = 0, vY=0 at £ =f and £ = |]. (14) 


We next expand W as a sine series of the form 
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v = >> C, sin ntt (15) 


thus satisfying the boundary conditions, ¥ = 0 at — = O and 1. With this expansion 


for V Eq. (12) becomes 
(D? — a’)W = (1 + af + BE’) D° C,, sin mat, (16) 
m=1 


and this equation is solved exactly to satisfy the remaining four boundary conditions 


on W. 
The general solution of Eq. (16) can be written in the form 


W = p oa {A cosh at + Bi” sinh at + A$”’E cosh at 
+ B,” sinh af + (n + ag + 6) sin mat 
4rm 
+ meta (a + 26¢) cos mt} (17) 


where A{”, B{”, A$” and B;” are constants of integrations and 


_ 48(a* — 5m*n") + (min’ + a’) 
~ (mx + a*)’ (18) 


The conditions that W = DW = 0 até = Oandé = 1 yield: 





yom — Amr 
mr +a 
(m) (m) __ — 86) 
aB;”’ +A, = (, a ma? + a? mr, 
A‘™ cosh a + Bi” sinh a + A$” cosh a + By” sinh a (19) 
— m+1 Amr 9 
_ ( 1) mn * + a (28 + a) 


A‘™asinh a + B}"’a cosh a + A” (cosh a + asinh a) 





8 
+ B,”’(sinh a + a cosh a) = (-9r"(a +B+at — :) ma. 


The solution of these equations is: 


= toma 

A; ee ee 3» 
mar +a 

- mr, , . ’ : : 

B;™ = faa,, + (sinh a + a cosh a)BZ — (sinh a)yZ}, 
A 
(20) 

ss — eee , —_— , ’ 

Ao” =- {(sinh’ aja, + (asinh a + a cosh a)B{, — (asinh a)yi,}, 


ae mr : . : 
BS™ = — {(sinh a cosh a — a)aZ + (a’ sinh a)BZ — (a cosh a — sinh a)yZ}, 
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where 


. 9 9 
A = sinh’ a —a, 


a 
—s mn” + a” s 
sa a Som a m+1 28) |. 9 
62 = a aw [ 1) (1 1+ + cosh a (21) 
, m+1 - ‘ 4a ° 
vi, =(-)D" "a+ 68+ a,) + —>— 3 (sinh a). 
mar +a 


With ¥ and W given by Eqs. (15) and (17), Eq. (13) becomes 
C,( 25 2 si rt - 7 a 14 0S b 
> nmr + a’) sin nr deta pn + a) i” cosh ag 
+ Bi” sinh aé + A}”’é cosh at + B,”’é sinh até 


4 ? 
+ (n + at + BE) sin mrt + eager — (a + 2B) cos we (22) 
mam + ¢ a” ) 


Multiplying this last equation by sin mzé and integrating over the range § = 0 tog = 1, 
we obtain 


>| os ad fn + (—1)"*' cosh aJA{” + [(—1)"*' sinh a]B,” 
\ 














= Ln’r’ ry a’ 
+ (—1)"| cos _ — ink | ed 
sh a mae a ssinha |A,; 
2a 
aa [ (- sinh a — I ;(1 + (—1)"*" cosh a) | = 
nner +a” 
1 a: a 
+ A2n.m + BYn.m + : le . (n*x”* + a’)’ 4a lk. = 0, (23) 
with 
K,, = C,.(m’n’* + a’)”’, (24) 
and 
z = 4 ee oe a aon if n=m; 
n,m 4 Yn.m = 6 An? x’ nx + a’ ’ — ? 
9 
Znm = 0 Yam = in| 2(n? = m’) bi mo ge z\, if —™ (25) 
oe i n + m even; 
a 4nm ! 1 2 | i ~ mM 
am” a ee. 2 2 2/2 eas - 2 i? 
n m Lx (n m°) mn +a + enelt. 


Equation (23) represents a system of linear homogeneous equations. It is the vanishing 
of the determinant of this system that determines 7 for any assigned a’. 
III. Numerical results. The equation determining 7 can be written as 
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| {It + (—1)"* cosh a]A” + ("sinh a 


|e +e 
n+1 2 b4 m 
+ (-1) [cosh a — Pts sinh a | : 
+ | (—9" sinh @ = es (1 + (— 1°" coh a) | =} 
nx + a’ 2 


| 
n = 1 2.2 2\3 5a. _ 

+ AXLn,m + a + 2 2 (n° + a) aT | — 0. (26) 
In solving the infinite determinantal equation (26), the determinant was terminated 
after three rows and columns. For a between —1.00 and —2.00 this seemed justified. 


For a < —2.00, the calculations indicate that a higher approximation, by terminating 





TABLE 1 
Cylinder ratios and values of a with their corresponding 8 and yp. 
Case R./R, d= R, = R, a B Bb 
I eh 1.345 1.035 75 2557 5058 
3.00 ‘ j . 7 . 
I 1.345 1.035 —1.0 .3410 .3410 
I 1.345 1.035 -1.5 .5115 .0115 
4.035 
II a 7 1.062 . 235 — .75 .0636 .3136 
II 1.062 235 —1.0 .0848 .0848 
II 1.062 235 —1.5 .1272 — .3727 
4.035 
III 355" 1.137 485 — .75 . 1276 .3775 
III 1.137 485 —1.0 . 1700 . 1700 
III 1.137 .485 —1.5 . 2551 — .2449 


the determinant after the fourth row and column, would be necessary. However, for the 
cases computed the convergence is fairly rapid as will be seen from Table 2. 

In the present calculations, three values of a were chosen which correspond to the 
experimental cases studied by Taylor; these cases are denoted by I, II and III in the 


TABLE 2 
Taylor numbers for assigned values of a, a and 8. 
Case a a B 7; T: T: T;/T: 
I 3.12 — .75 . 2557 2464.84 2462.38 2453.04 . 9962 
II 3.12 — .75 .0636 2668 . 86 2662.94 2652.28 . 9960 
III 3.12 — .75 .1276 2597 . 24 2592.75 2582.54 .9961 
I 3.12 —1.0 .3410 2885.15 2878.25 2867 . 97 . 9964 
II 3.12 —1.0 .0848 3276.10 3256.78 3244.32 . 9962 
Ill 3.12 —1.0 .1700 3134.83 3120.86 3109.20 . 9963 
I 3.20 —1.5 .5115 4383 .32 4328.23 4316.76 .9973 
II 3.20 —-1.5 .1272 6018.53 5761.04 5747.24 .9976 
III 3.20 —-1.5 2115 5353.70 5198.56 5181.98 . 9968 
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subsequent work. For these three values of a, Eq. (8) leads to three values of 6 and u. 
Case II gives the smallest separation of the cylinders and hence should be approximated 
best by the theory. Case I is for the widest separation of the cylinders. 

Table | gives the ratio for each case, the value of d and the three values of a with 
their corresponding 8 and uz. 

Table 2 lists the Taylor numbers calculated for the nine values of 8. Here the sub- 
script on 7 indicates the order of the determinant which was made equal to zero in the 
computation of 7’, . The final column gives the ratio of 7 to T, and provides a measure 
of the convergence of the method. 

The critical values of a for these cases agree exactly with those used in the first 
approximation of y/d. Therefore it was possible to use some of Miss Elbert’s com- 
putations. This was accomplished by writing the constants of integration in the forms 


Ao = AD’. 
(m) m mrp f m+1 8 
By” = By’ + —j\(a + (-1) cosh a) —33->—— 
A mrmr-+a 
a , nok 2 . 4(a® — 5m*r’) Pome 
+ (a — (—1)”*' sinh a) ——355 xe — (-1)""'sinh a, 
(mr + a) ) 
- mrB},..9 ; ee 8 
so” = AP? — —4(sinh’ a + (—1)"*'a’ cosh a) —3-3——3 
a mae +a iaied 
(27) 
wwe — : 4(a®” — 5m’) a ; 
+ (sinh? a — (—1)”*'asinh a) ——35 sa- — (—1)"*'asinh a?, 
(mr + a) ) 
rer a mrB),. ae 
B” = By + A ‘(sinh a cosh a — a + (—1)""' sinh a 


8 1a’ — 5m*n’) 
m+1 
— (—1)”’ a cosh a)( = 33 2 
mr +a (mr +a) J 
\ 


ee ee 8 ve : \ 
+ (—1)”* a’ sinh a — > :) — (—1)”"" (a cosh a — sinh a)?, 
mar+a J 

where the subscript f denotes the constants from the first order calculation. Obviously 
the corresponding value of u varies. When a has the value —2.00 or smaller the critical 
value of a no longer agrees with that found in the first approximation but rather appears 
to be consistently smaller. 

A comparison of the values derived in this paper with those obtained by both Taylor 
and Chandrasekhar is made in Table 3. The experimental results are given in units 
of 2,/v with possible deviations of one or two per cent in the range considered. Equation 


(28) was used to convert the computed Taylor numbers to the unit of measure Q,/7; 
(1) _ —70 — R,/Ri) (28) 
v id*(1 — uR2/R)) 


In this table the first column lists the case; the second, the value of yu; and the last four 
columns the values of 2,/v from Chandrasekhar’s first order calculation, from the 
nearest second order calculation, and from Taylor’s theoretical work and lastly from 
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TABLE 3 


Comparison of values of 2;/v derived from the first and second order calculations using Chandrasek- 
har’s method, from Taylor’s theoretical work and from Taylor’s experiments. 


Case rv Q1/veirst 21/Vsecond 2; /vpaylor 21/vexperiment 

I . 5058 68.6 71.3 77.4 74.2 

I . 3410 34.2 36.3 36.6 39.4 

I .0115 27.3 31.0 27.8 30.5 

II .3136 204.9 207.5 207 .2 203.0 

II 0848 190.6 194.0 193.9 193.8 

II — .3727 197.2 206.0 203.7 212.2 
III .3775 79.9 81.7 81.5 86.8 
III . 1700 70.2 72.6 72.8 74.6 
III — .2449 66.8 72.2 71.6 73.1 


Taylor’s experiments. Both the Taylor values and the first order values had to be inter- 
polated. 

It is evident from Table 3 that in all cases the second order calculation gives a 
value that is larger and closer to the experimental value than the first order calculation. 
It is also true that for the particular cases considered the second order value is more 
often nearer to the experimental values than is Taylor’s value. As a + — © the Taylor 
theoretical values will deviate more and more from the experimental values. It seems 
that a third order calculation would add a correction of the right magnitude. However, 
this would probably require the inclusion of the second order terms in p. 

The author is indebted to Dr. 8S. Chandrasekhar for suggesting this problem and 
for his valuable advice. Thanks are also extended to Dr. Chandrasekhar and Miss 
Donna Elbert for the use of their first order computations. 
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BOOK REVIEWS 


Linearized theory of steady high-speed flow. By G. N. Ward. Cambridge at the University 
Press, 1955. xv + 243 pp. $6.00. 


This excellent monograph is a revised version of an essay submitted in competition for the Adams 
Prize for 1949-50 in the University of Cambridge. The published version has been brought up to date 
so that it contains most of the important work in the field up to, and including 1953. The book is divided 
into three parts; in the first part the linearized equations of motion are derived, general solutions of these 
equations for subsonic and supersonic flow are presented, and the boundary conditions, aerodynamic 
forces, as well as uniqueness and flow reversal theorems are discussed. The second part deals with special 
methods of solution for two-dimensional and three-dimensional wing problems for both subsonic and 
supersonic flow, although much more emphasis is placed on the latter case. This portion contains sections 
on subsonic flow past thin bodies; supersonic flow past nearly plane wings, with full use made of the 
more sophisticated techniques such as “finite part of an integral”; conical fields in supersonic flow where 
the methods of complex function theory are utilized; and finally a section on the application of operational 
methods to problems of supersonic flow. The third and final part deals with slender-body theory, a field 
with which the author’s name is most intimately connected. The monograph concludes with an inclusive 
bibliography of nearly 200 references. 

Vectorial notation has been used extensively throughout and proves to be very useful in many 
instances, although this reviewer cannot wholly agree that it results in a more satisfactory approach 
than to work in terms of potentials. The reviewer agrees with the author that a knowledge of the elements 
of compressible fluid theory is required of the reader. Certainly, the book recommends itself highly for 
applied mathematicians, theoretical aerodynamicists, and research workers engaged in the study of prob- 
lems of high speed flow. Moreover, there is no reason why it should not also make an excellent text for 
a graduate course in linearized compressible flow theory. One can only conclude by pointing out that 
this monograph continues to maintain the high standards in content, as well as in format and printing, 
which have already been established in this series of Cambridge Monographs in Mechanics and Applied 


Mathematics. 
RonaLp F. ProspsteIn 


Mathematical foundations of quantum mechanics. By John von Neumann. Translated from 
the German edition by Robert T. Beyer. Princeton University Press, New Jersey, 
1955. xii + 445 pp. $6.00. 

This is a translation of von Neumann’s well-known book, Mathematische Grundlagen der Quanten- 
mechanik which first appeared in 1932. The translator’s preface is dated 1949, the long interval which has 
elapsed before publication not being explained. Since the original edition appeared, there have been 
important advances in Quantum Mechanics, notably in connection with the quantization of wave fields. 
No attempt has been made to bring the book up to date in this regard, and the translation sticks very 
closely to the original German edition. Even the references to text-books have not been brought up to 
date. 

The translation is somewhat too literal in many places, and the translator often employs phrases 
and words which, if not absolutely incorrect, have at least a strange and somewhat “foreign”’ flavor to 
them. Nevertheless, the meaning is always clear. One slight difference with the original which most people 
will probably consider an improvement, is that the footnotes have been placed in the body of the text 
instead of at the end of the book. 

The translator and publisher have performed a service in making this classic available to a wider 
circle of English-speaking readers. It remains indispensable to those who desire a rigorous presentation 
of the foundations of the subject. The price seems rather high in view of the quality of the printing and 


the paper binding. 
A. F. STEVENSON 
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ON UNIQUENESS IN THE THEORY OF PLASTICITY* 


BY 
D. C. DRUCKER 
Brown University 


Summary. The fundamental definitions of work-hardening and perfect plasticity 
have far reaching implications with respect to uniqueness of solution for elastic-plastic 
bodies. Satisfaction of the basic postulate, that in a cycle work cannot be extracted 
from the material and the system of forces acting upon it, guarantees an existing solution 
to be stable but not necessarily unique. Uniqueness follows for the usual linear relation 
between the increments or rates of stress and strain and also for combinations of such 
linear forms. Conversely, lack of uniqueness results for an elastic-perfectly plastic body 
when, for example, the maximum shearing stress criterion of yield is employed with the 
Mises flow rule. 

Uniqueness [1]{2]{3][4]'. Plastic stress-strain relations for work-hardening materials 
are strongly path dependent. It is, therefore, not reasonable to expect that a given set 
of final boundary values will give a unique solution for stress and strain in the interior 
of a body. A complete specification of the history of the applied surface tractions, body 
forces, and surface displacements will generally be required. It is simpler then to speak 
of a body under existing surface tractions 7; , body forces F, , displacements u, , stresses 
o;; , and strains e;; . The question is then whether the stress and strain rates, o{; and €/; 
are determined uniquely by the rates of change of the applied forces and displacements 
T: , Fi, and uj. 

The terminology of plasticity is developing continually and is much a matter of 
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Fig. 1. Tresca criterion and associated Fig. 2. Mises criterion and associated 
flow rule. flow rule. 








*Received Feb. 23, 1955. The results presented in this paper were obtained in the course of research 
sponsored by the Office of Naval Research under Contract N7onr-35801 with Brown University. This 
note is a more formal presentation of a talk given at the Applied Mathematics Seminar of the University 
of London November 1954. 

1Numbers in brackets refer to the bibliography at the end of the paper. 
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Fig. 3. Convexity and normality. 


individual preference. A brief description therefore will be given of some of the terms 
and concepts employed here. The state of stress o;; at a material point of the body 
may be represented by a point o,; plotted in a stress space. Coordinates of the point 
are the components of stress. The vector from the origin to the stress point is called 
the stress vector. If, for example, the normal stresses o, and o, are the only non-zero 
components of stress, the familiar two-dimensiona! plots of Figs. 1 and 2 are useful. 
In general there are nine components of stress, six of which are independent because 
cross shears are equal o,; = o;; . The general stress space is nine-dimensional and is 
shown symbolically as in Fig. 3. 


“ey ~ "eage 
pe, eae. illic ek 
a b . 
4a 4b 


4c 4d 
Permissible : b—-a for “e'? 
poth ob b—-a — a—eb for "e'? 


Fig. 4. Permissible paths. 


When the stress at a material point changes, the stress point moves about in stress 
space. The region in which only elastic changes in strain occur is bounded by a surface 
called the yield surface. Figures 1 and 2 show the two most commonly assumed yield 
criteria. For a work-hardening material, a change in stress which causes the stress 
point to move outside the existing yield surface is termed loading. Both plastic and elastic 
deformation will be produced and a new yield or loading surface will be established. 
Subsequent yield surfaces need not resemble the original in appearance. Again, changes 
in stress which move the stress point about inside the region bounded by the new surface 
will cause elastic changes only in strain. Perfect plasticity is a limiting case of work- 
hardening. The yield surface is fixed in stress space as the material cannot carry stresses ) 
“above” yield. The stress point may not move outside the surface and plastic deforma- 
tion will occur for points on the yield surface. 
If strain coordinates are superposed on the corresponding stress coordinates, the 
strain rate may be exhibited as a free vector on the same diagrams. It is generally de- 
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sirable to place the strain rate vector ¢/; at the stress point o;; to emphasize that the 
strain rate is associated with both o;; and o{,; . For convenience and simplicity the plastic 
component of the strain rate, ¢«/? , alone will be plotted. The scalar product of the stress 
vector and the plastic strain rate vector o, ;€{; is rate of plastic work or energy dissipation. 

Consideration of a homogeneous state of stress and strain in a plane and the Tresca 
yield criterion with its associated flow rule, Fig. 1, illustrates the well known fact that 
too much should not be expected in the way of uniqueness. At point A, many plastic 
strain rates ¢’” are associated with a given state of stress [5]. On side AB, many states 
of stress are associated with a given plastic strain rate. Any such lack of uniqueness in 
the small will show up in the study of uniqueness in general. 

Uniqueness proofs ordinarily follow a standard pattern [1]-[4] which will be reviewed 
briefly before introducing a new point of view. Two solutions a and b are assumed; 
°o!, , “ef; and °o!; , °ef; corresponding to 7; on the boundary A, , to u{ on the remaining 
boundary A, , and to F{ in the volume v. The theorem of virtual work is then employed 


(repeated indices denote summation): 


[ Tu, dA + [ Fru. in [ ote. mt (1) 
A ve 


ae, vo 


The starred quantities are related through equilibrium and the unstarred are compatible. 
There need be no relation between the two sets of quantities. The difference between 
the two assumed states a and b, therefore, can be substituted in Eq. (1) although 
*o’, — °o!; may not produce “*e/; — °e/; . Substitution gives 


Oi; 
0 = il [°o?; — °of;)[°el; — °ef;] dv (2) 


because °7” = °T’ on A7, “us = “ui on A, and °F = °F!. 

If it can be shown that the integrand in (2) is positive definite, uniqueness is proved 
to the extent at least of either o/,; or e{; having but one possible value at each point of the 
body. In what follows, the term uniqueness will be used without qualification. As a 
first step, the strain rates are ordinarily decomposed into their elastic and plastic portions 


ef, = 5 + ef} (3) 


ii ii 
and the integrand is written 
[°o!; — Pols )[*el} — “elf] + [ots — “ots )[Pet? — *ef7). (4) 
The first term is positive definite for both linear and nonlinear elasticity. If then the 
second term is positive or zero, uniqueness is established. For a rigid-plastic material 
the first term is identically zero and this sufficiency proof requires the second term to 
be positive when °o/; ¥ °ot,; or “el; # "Gs 
At this stage it seems appropriate to introduce the alternative point of view which 
avoids all mathematics for a wide class of stress-strain relations. 
Mechanical-thermodynamic postulate. The fundamental definition of a work- 
hardening material which has been advanced previously [5][6] is that over a cycle no 
work can be extracted from the material and the system of forces acting upon it. An 
alternative and more precise set of statements refers to an external agency which applies 
and removes a set of forces to the already loaded body. The external agency must do 
positive work during the application of force. Over the cycle of application and removal 
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of force the work done by the external agency must be positive if plastic deformation 
occurs and will be zero if elastic changes only take place. The definition covers the usual 
theories and any theory of plasticity based upon a work-hardening mechanism. Slip 
theory of Batdorf and Budiansky [7] which is based upon a work-hardening shear 
mechanism, or polycrystalline theory based upon single crystals which individually 
satisfy the work-hardening requirement [8] are good examples. 

The definition for an ideally or perfectly plastic body is similar but the work done 
by the external agency may be zero when plastic deformation takes place. 

It follows from consideration of a homogeneous state of stress that the yield and 
subsequent loading surfaces are convex and that the plastic strain rate vector is normal 
to such surfaces at a smooth point or lies between adjacent normals at a corner [5], 
Figs. 1, 2. These results may be stated in mathematical form, see Fig. 3, as 


(o;; — o*)e? > 0, (5) 


oie; 2 0 (work-hardening), 
= 0 (perfectly plastic). (6) 


The equality sign in (6) applies to a work-hardening material only for the trivial case 
of e{; = 0. 

Permissible stress paths and a uniqueness proof. Satisfaction of the basic thermo- 
dynamic postulate leads inevitably to convexity of the yield and subsequent loading 
surfaces and also to generalized normality of «/? . These results must be obtained no 
matter how elaborate or how simple the subsequent physical reasoning and mathematical 
development. The postulate might seem also to be powerful enough to insure uniqueness. 

According to the concept of work-hardening or of elastic-perfectly plastic bodies, 
the external agency must do work to change one stress state in a body to another. In 
other words, constant external forces cannot, in themselves, provide the necessary 
addition to the sum of the stored and dissipated energy. This does give a certain type 
of uniqueness. If a state of stress is attained at a given system of loads, it is stable. 
Solutions for elastic-perfectly plastic bodies are therefore unique. However, uniqueness 
of solution for work-hardening bodies is more general. Stability does not necessarily 
rule out the possibility of two different states being reached as the external loads are 
changed by very small increments. There need not be any permissible path between 
the two states corresponding to no change in load. In a sense an insurmountable local 
energy barrier may exist. 

If, however, the system is truly linear in the increments or rates of force, displacement, 
stress and strain a permissible path does exist. A combination of two alternative solu- 
tions would also be a solution. The stability and the uniqueness questions are the same. 
Thus the fundamental postulate insures uniqueness for linear and non-linear elasticity. 

The incremental linearity of the usual plasticity theory and the extended incremental 
linearity obtained by combinations of linear forms [4][10] are obviously not linearity 
in the true sense. Plastic action is an irreversible process. Coefficients appearing in the 
linear expressions relating the rates of stress to the rates of strain depend upon whether 
the stress changes produce elastic-plastic strain or purely elastic changes. There will be 
variation from point to point of the body for this reason as well as for the more usual 
ones applicable to the plastic strain rates. Similarly the differences in strain rates be- 
tween the two assumed solutions a and b will be related linearly to the differences in 
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the stress rates. The coefficients will vary from point to point in the body depending 
now upon both solutions. 

A pseudo-material may be defined with the same linear incremental stress-strain 
relations which apply to the difference between a and b but with the added physically 
improper assumption of reversibility. The question of uniqueness for the real material 
is then replaced by a question of stability of solution for the pseudo-material. 

The work done by the external agency in going from a to b is the same as for b toa 
for this pseudo-material. A permissible path from a to b or b to a always does exist, in an 
extended sense, in the actual material Fig. 4. The work done by the external agency, 
therefore, is positive at each point. Two solutions cannot exist for the pseudo-material 
and uniqueness follows for the real material. It should always be kept in mind that 
uniqueness may be limited by conditions as in Fig. 1 and that rates of change of overall 
geometry are assumed negligible. 

The preceding paragraphs are a proof, in words, that the fundamental postulate 
does insure uniqueness for the linear incremental theories of plasticity which enjoy 
such popularity and for the more recent type of theory which employs combinations 
of incrementally linear mechanisms [4][7][8}[{9][10]. 

A return to the more conventional uniqueness proof will, perhaps, help to clarify 
the considerations of irreversibility and permissible path. An infinitesimal time may 
be taken to elapse for each of the two assumed stress rates separately. The state of 
stress at a point in the body will change from the existing stress o;; to the stress points 
a or b depending upon which assumed alternative solution is followed. Some typical 
examples are drawn in Fig. 4. If a and b are both elastic changes, Fig. 4a, expression (4) 
is positive. If b is elastic and a is elastic-plastic, Fig. 4b, the second term of expression 
(4) is positive from (5) and (6) and the entire expression is positive. Reversing a and b 
does not alter this result. When a and b both represent elastic-plastic changes, Fig. 4c, 
and the incremental relation is linear 


ey - Ai jth ’ (7) 


the difference between a and b will satisfy (6) and (4) will be positive definite. The 
coefficients H;;,, are functions of stress and may depend upon the strain and the history 
of loading [12]. 

It is instructive here to look at the difference between states a and b from the point 
of view implied by the fundamental definition of work-hardening. Call the infinitesimal 
time unity so that the two assumed states of stress are o;; + °o/; and o;; + °o/; . The 
corresponding strains are €;; + “e/; and e;; + °e/; . It may be possible to go from strain 
state b to strain state a by going in a straight line in space from b to a. The step 


*o’, — °of; may be considered then as applied by an external agency which produces 
the corresponding strain “e{; — °e{; . The fundamental postulate then gives 
b b 
[°ot; — “ois Je; — “eis] > 0 (8) 


and uniqueness is established. Figures 4a and 4b are evidently in this category as is 
4c if the linearity relation (7) applies. Note, however, that all paths are not permissible. 
In Fig. 4b, for example, the path must go from the elastic to the elastic-plastic domain. 
It will generally be true that for some regions of the body the permissible path will be 


from b to a and for other regions from a to b. 
If a corner exists as at A of Fig: 1, Fig. 4d, and the assumption is made that two or 
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more independent linear mechanisms operate, uniqueness again follows in an obvious 
way by treating separately the strain components each mechanism produces. Here 
the path may be from a to b for one component and b to a for another. 

As can be seen, the permissible paths are somewhat peculiar and linearity has a 
rather unusual meaning in incremental plasticity theory. The pseudo-material which is 
linear and reversible has no physical validity. Nevertheless, satisfaction of the funda- 
mental work-hardening requirement plus incremental linearity in the ordinary or ex- 
tended sense insures the development of a stress-strain relation which will give unique 
solutions. 

If the stress-strain relation is incrementally non-linear [5], the fundamental definition 
assures stability of solution. However, inequality (8) does not necessarily hold because 
there is no stress path from b to a in the neighborhood of ¢;; which will take the strain 
state b to the strain state a or vice versa. 

The absence in plasticity of the overall uniqueness which holds in elasticity is closely 
connected with the impossibility of such paths for the real or for a reversible pseudo- 
material. Given any finite change in 7’; , for example from zero, it might be thought 
uniqueness could be shown as in elasticity through an inequality like (8) in which the 
primes are dropped. The point is that in general, and there are exceptions, it is not 
possible to change strain state b to state a by going from stress state b to a along a 
path for which the virtual work term [*e;; — °o;;][*e:;; — °e:;] has the same sign as the 
actual positive work done by the external agency. 

Perfectly plastic bodies [11]. As has been stated, the problem of stability of solution 
and of uniqueness are essentially the same for elastic-perfectly plastic bodies. The 
definition of perfect plasticity does therefore, assure uniqueness. A conventional type 
of proof is sketched below for comparison. 

The yield surface for ideally or perfectly plastic material is fixed in stress space. 
The state of stress cannot be outside the yield surface and plastic deformation occurs only 
when the stress point is on the surface. Normality holds at a smooth point and extended 
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normality applies at a corner. Considering a corner as representing the sum of independent 
mechanisms, the normality condition applies to each independent yield curve and there 
need be no distinction in the treatment of the two cases. Although the plastic term in 
expression (4) may be zero, the elastic is always positive unless *o/; = °o/; . Uniqueness, 
in this sense, is established for the elastic-plastic case but not for the rigid-plastic where 
the elastic term is identically zero at all times. 

Tresca yield criterion and Mises flow rule. It has sometimes seemed convenient to 
adopt a maximum shearing stress criterion of yield and to couple it with the Mises 
flow rule, Fig. 5, in solving problems in perfect plasticity. The deviation of ¢{} from the 
normal to the yield surface is the significant feature. It means that in a cycle work can 
be extracted from the material and the system of forces acting upon it. This violation 
of the fundamental postulate would seem to imply a lack of uniqueness which in turn 
casts some doubt on the validity of solutions obtained by such starting assumptions. 
In terms of the usual uniqueness proof, expression (4) can be made negative. Although 
the fundamental postulate or equivalently the positive definiteness of (4) provide suffi- 
cient conditions for uniqueness, their necessity does not follow. It might be true for all 
problems that a negative value of (4) or of work done by the external agency in one 
region of the body would be more than balanced by positive values elsewhere. 

Proofs of necessity are difficult and often tied in with questions of existence of solu- 
tions. Such a proof will not be given. The following demonstration by means of a par- 
ticular solution will serve to show lack of uniqueness for the material of Fig. 5 and make 
necessity quite plausible. 

Consider a tube with a very thin wall subjected to interior pressure p and acted 
upon by an axial elastic spring under tension F and by a compressive force 7 which 
relieve some of the axial stress induced by the interior pressure, Fig. 6. The stress point 
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Fic. 6. Model to illustrate lack of uniqueness. 


will be taken as point S, Fig. 5. Instability or lack of uniqueness shows up immediately. 
Consider a compressive surface traction rate 7’. Equilibrium, compatibility and stress- 
strain relations are satisfied by a purely elastic decrease in height. Both the spring 
force and the axial tension in the tube wall would then decrease. If plastic deformation 
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is assumed to occur, the situation is quite different. The plastic strain in the tube has 
a component of axial compression. This shortening of the tube can decrease the force 
in the spring to any desired extent and so increase the axial tension in the tube. It is 
easily seen that the lack of uniqueness is equivalent to an instability of solution. At 
T = 0, as plastic deformation takes place, the force in the spring continues to drop as 
the axial tension builds up in the tube until the strain rate vector has no axial component. 

This lack of uniqueness does not occur if the fundamental definition of perfect 
plasticity is adhered to. The flow rule appropriate to the Tresca yield criterion gives a 
normal plastic strain rate vector which has no horizontal component at S. Work cannot 
be extracted; o/,;¢/? cannot be negative as in Fig. 5; 7’ produces the unique result that 
F and the axial stress in the tube decrease elastically. 

It might well be argued that the example of Fig. 6 is not a good one because the change 
in tube diameter, which has been ignored, will increase the stress and there is an insta- 
bility in this sense. However, in principle at least, the stresses could be set up in a flat 
sheet to obviate this trouble. It is the possibility of work extraction represented by the 
angle between the plastic strain rate vector and the normal to the yield surface which 
is responsible for the instability or lack of uniqueness. 
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ON THE VIBRATION OF ELASTIC BODIES 
HAVING TIME-DEPENDENT BOUNDARY CONDITIONS* 


BY 
J. G. BERRY AND P. M. NAGHDI 
University of Michigan 


1. Introduction. In the classical treatment of the small vibrations of elastic bodies, 
the displacements are assumed to have the form’ 


U.(z, t) = uii(x)qn(t), (x = 2 ’ Z2 ’ Zs), (1) 


where the u?(x) are the so-called normal modes and the q,(¢) are the generalized co- 
ordinates. The boundary conditions, which are homogeneous, serve to determine the 
mode shapes and the secular equation governing the natural frequencies. If, however, 
the boundary conditions are non-homogeneous, the classical method usually fails since, 
in general, the boundary conditions and the displacements involve different functions 
of time.” 

Among the various methods which may be employed in the study of the motion of 
a vibrating elastic body having time dependent boundary conditions, the integral 
transform techniques and the use of Lagrange’s equations,’ for certain problems, are 
perhaps most familiar. An alternative approach is to transform the original problem 
into an equivalent forced vibration problem with homogeneous boundary conditions; 
this approach is adopted here. The basic idea involving the transformation of variables 
so as to remove the non-homogeneous boundary conditions, is well known [1l, p. 277]. 
However, this procedure, in connection with vibration of elastic bodies having time 
dependent boundary conditions, was first used by Mindlin and Goodman [2] in dis- 
cussing the motion of beams and later by Herrmann [3] for problems of rods. 

The present paper contains a method of solution for the motion of any finite, vibrating, 
elastic body having arbitrary time-dependent boundary conditions and initial con- 
ditions. The method used involves a transformation of the type employed by Mindlin 
and Goodman [2]. However, the resulting forced vibration problem is discussed without 
reference to stress strain relations and therefore the procedure is quite different from 
that contained in [2] or [3]. Since the stress strain relations are not explicitly specified 
the results are readily applicable to all linear theories of elasticity including the special 
theories of rods, beams, plates and shells; in each case, it is only necessary to introduce 
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1Throughout this paper, the indices (subscripts or superscripts) 7 and j take the values 1, 2, and 3 
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upon introduction. The indices i and j are reserved for quantities which have tensor character while 
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over all the values the indices may take. Whenever one of the repeated indices is placed in parenthesis, 
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2An important exception occurs when the boundary conditions are simple harmonic functions of 
time and only the steady state solution is sought. In this case, no difficulties are encountered. 

*Some remarks concerning the use of Lagrange’s equations are given in Sec. 5. 
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the appropriate stress strain relations. This, together with some of the features of the 
procedure, is illustrated by means of an example in Sec. 4. 

2. Statement of the problem. Let V and S be the volume and the boundary surface 
of a region R (not necessarily simply connected) occupied by an isotropic, elastic medium.* 
Further, let U; , o;; , X; and p denote, respectively, the displacement vector, the stress 
tensor, the body forces, and the mass density all referred to a set of rectangular cartesian 
coordinates x; . Then, the equations of motion for the medium are 

Oi3.5 + Xi = p OU, /at, (2) 
where comma denotes partial differentiation with respect to the space variables. 

Let B, denote a system of spatial differential operators such that when they are 
applied to appropriate components of the displacement vector (written as B,[U, , U2, 
U,]), these quantities take on prescribed values on the boundary S and will be assumed 
to have the product form g(x) f(t). Thus, if there are p boundary conditions, r(r < p) 
of which are non-zero, the boundary conditions may be written as 


B,[U, ; U2 , Us] = g(a) f'(O8e on S, (3) 
where 6; is the Kronecker delta, and the indices k and / have the range k = 1, 2,3 --- p, 
l= 1,2,3--- r. It is important to note that the notation B,[U, , U., U3] does not imply 


that the operators B, are applied to the displacement vector U; ; the B, are applied 
only to appropriate components of U; [see for example Eqs. (23)]. 

We now wish to determine a set of displacement functions U; which satisfy the 
equations of motion throughout R, the boundary conditions on S, and the initial con- 
ditions in RF for allt < 0. 

3. Development of the method. We introduce the transformation 


Ua, ) = v(x) f'() + é(2, O, (4) 


where v; denotes the /th displacement vector which will be defined presently. If we now 
. A . l l l el — C 

require that the functions v; take the values B,[v; , v2 , v3] = ged) on S, it then follows 

from (3) and (4) that the displacements £; must satisfy the homogeneous boundary 


$ 


f,, £,&] = OonS. 


; 


r 
| 


conditions B,| 

A suitable choice of v; may be made in various ways. Mindlin and Goodman [2] 
represented v; as polynomials with arbitrary constants which were then adjusted to 
give the desired results. A more general method, which reduces to that of [2] for certain 
classes of problems, is to choose v; to be the displacements which satisfy the following 
system of statical problems: 


= 0) in Rf, 


“) 


(5) 


l l l \ecl Y 
B,[v; , v2 , 3] = ge(2) be) on S, 


where 7!; is that portion of ¢;; which is obtained from the displacements v; . Equations 
(5) represent a system of (l) separate problems, each of which has been obtained from 
the original problem [Eqs. (2) and (3)] by setting X; = 0, 0°U;/at’ = 0, replacing o;; by 
ri; , and setting all f(t) = 0 except that f‘(t) = 1 for each value of / in turn. Thus, each 


4Actually, the medium need not be assumed to be isotropic. In fact, the subsequent analysis is valid 
whenever the medium possesses a strain energy function which is positive definite. 
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of the r sets of v; satisfy the same equilibrium equations but different boundary con- 
ditions. In the special case where the boundary conditions are all displacement conditions 
(no surface tractions), this procedure coincides with that used by Herrmann [3]. 

Due to the manner in which the boundary conditions appear in (5), it is possible that 
the region F is not in statical equilibrium. For example, consider R to be a beam pinned 
at one end only and acted upon by an oscillating couple such that the beam rotates 
about the pin. The corresponding statical problem characterized by (5) is a beam pinned 
at one end only and acted upon by a constant couple; it is evident that the beam is not in 
statical equilibrium. This difficulty is easily resolved if we interpret statical equilibrium 
according to D’Alembert’s principle. Thus, if the region R is undergoing rigid body 
motions due to the tractions on S, it is sufficient to add to the first of (5) those body 
forces necessary to insure equilibrium. We shall denote these body forces by Q)(x). In 
general, there will be one set for each of the (/) separate problems. The necessity for 
introducing the forces Q! is really a consequence of the manner in which the v; were 
chosen. 

In what follows, it will be assumed that the displacements v; are known. In view 
of the character of the transformation (4), the displacements £; are the solutions of a 
forced vibration problem with homogeneous boundary conditions. In keeping with the 
classical approach to such problems, we assume £; to have the form 


E(x, t) = ui(x)q,(2), (6) 


where wu; are the normal modes of the associated free vibration problem and the functions 
q,(t) are the generalized coordinates, as yet undetermined. 

Before proceeding further, it is advantageous to deduce certain relations associated 
with the normal modes. To this end, let 7;(z, ¢) in the form 


mn = > ui(z) exp [(—1)'any] (7) 


n 
be the displacements which satisfy 


Si; ee d°n;/dt? in R, 
(8) 


B,{m » Ne » ns] = () on S. 


Equations (8), which serve to determine the normal modes and the natural frequencies 
w, , are obtained from (2) and (3) by putting X; = 0, replacing U; by 7; , o;; by S,; , 

and setting all f'(t) = 0. 

The kinetic energy associated with 7; is given by 
T = 3 | o(6n./ad(n./ab) aV (9a) 
- 
and in view of (7), the maximum kinetic energy corresponding to the nth normal mode is 
T! = het i] ouu dV. (9b) 
- 


The normal modes, in addition to (9b), also satisfy Clebsch’s [4] orthogonality conditions 


[ puyu; dV = 0, (n ¥ m). (10) 


“Vv 
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With the aid of (7), the first of (8) can be written as 

D St. = —pwrul, (11a) 

where S;; is that portion of the stress tensor o;; which is determined from the nth normal 


mode. If Eqs. (11a) are now multiplied by wu”; and integrated throughout R, there results 
for each n 


[. S,,dV = —wy | putu™ dV (11b) 
V V 


which, by (9b) and (10), becomes 


I 
o 


/ u"S,, dV forn * m 
P (12) 


= —27’ forn = m. 


As will be seen shortly, Eqs. (12) play an important role in the determination of the 


generalized coordinates. 
Substituting (4) into (2) and with the aid of (5), (6), and (8), Eq. (2) will appear as 


tit (O + Stat) + X; = pr d’f'/dl + pul a’q,/at’ , (13) 


where, according to (5), the first term on the left hand side of (13) is either zero or 
— Qif'(t), depending on whether the region R is in statical equilibrium or not. The 
differential equation governing q,(¢) is now obtained by multiplying (13) by wv’, integrating 
over the volume V of R, and using (12) with the result 


O°”, 


2 
——- > = 3 of i r her! Ne , 
ot’ + Win) In [ ul X. — 0 i> Qif' Ja / [io d}\ (14) 


the general solution of which is 
t 
g(t) = A, SiN Wat + B, COS wey) t + + | P,(7) sin w(t — 7) dr. (15) 
(n) 0 


In (15) A, and B, are constants and P,(¢) represents the right hand side of (14). If the 
initial conditions are 


Uz, 0) = U,(2), dU,(z, 0)/at = U{,(2), (16) 


then, with the aid of (10), the following expressions for A, and B, are readily deduced: 


A, [ mal i= ((2LO |] av / oe [ ou u” dV, 


2. = / pur[U so -vsolav / uu” dV. 
Vv 


This completes the formal solution of the problem. By way of recapitulation, the steps 
leading to the solution of the problem under consideration are: (a) assume that the 


(17) 
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displacements are given by (4) and (6); (b) determine the displacements v! by solving 
the system of statical problems defined by (5); (c) solve the free vibration problem 
obtained by setting X¥; = 0 and f‘(#) = 0 in (2) and (3); and then (d) evaluate the 
integrals appearing in (15) and (17). 

It may be noted that whenever the boundary conditions involve simple harmonic 
functions of time, the above procedure can be modified slightly. In such cases it is con- 
venient to choose v} to be a set of steady state solutions corresponding to each time 
dependent boundary condition. This modification simplifies the results somewhat. 

4. Anexample. Consider a rectangular beam of length L, thickness h, and width b. 
Let the beam be referred to a set of cartesian axes x; with the length measured along 
x, , while x, is positive downward; the origin of x; is placed at the centroid of the left 
end section. The beam is assumed to be simply supported and is driven by a step- 
moment at the left end; the step-moment is the product of a constant moment M, and a 
unit step function in time. 

The equations which govern the motion of such a beam are obtained by assuming 
that the non-vanishing stresses are o,, and o,. and that these are functions of z, , x, and 
t only. We also assume that the bending stress ¢,, varies linearly across the thickness. 
A set of displacements which are consistent with the last assumption, are 


U, sad T2V(Xi ’ t), U, - w(x, ’ t). (18) 


The equilibrium equations and the stress displacement relations which may be readily 
obtained from the general equations of elasticity (see, for example, Ref. [5]) are 


aM ay 
Oz, as Q — pl at’ ’ 
(19) 
aQ _ 4 dw 
ae, 4 Pe 
and 
_ py ov eit ( 2) 
M = EI aft Q = pkAl y+ az,)? (20) 
where the bending moment M and the shear force Q are defined by 
+h/2 +h/2 
M = > | 01122 dx, ’ Q = b O12 dx, . (21) 
—h/2 —h/2 


In (19) and (20), E is Young’s modulus, u is the shear modulus, J = bh*/12, A is the 
cross-sectional area, and k is a constant whose physical significance has recently been 
discussed by Mindlin and Deresiewicz [6], but need not concern us here. The substitution 
of (20) into (19) yields the following two differential equations for y and w. 


ay :A( ou) ay 
EI 573 — ukAly + 5°) = ol 5, 


ox; 


(22) 


ov 4 He) 44H 
ura(2e ¥ ozi/ pA 
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which are generally attributed to Timoshenko [7]. Also for the example stated, the 
boundary conditions are 


B,[z.¥, w] = w0, ) = 0, Balte¥, w) = EL wee = M,f(t), 
O41 





dV(L, _ 


B;[z.¥, w] = w(L, t) = 0, B,|z.¥, w]) = El i 0; (23) 
(0 fort <0 
f = 4 
ll fort> 0. 
According to Eqs. (4) and (6) we take y and w as 
rey =v, f() + uigq,(d, 
(24) 
w = vf(t) + uq,(0, 
where »; and v, satisfy the equilibrium equations 
d’(xz"'v,) ee dv, | 
GI 3 — wh As (23'0,) 2) = 0, 
EI dr? 7 Ay (aa M1 ams ( 
(25) 
U(xz'v “Ve 
Aaa) 5 So 0, 
dx, dx; 
and the boundary conditions 
17 Uxz 'v,(( 
B,[v, , v2] = v.(0) = 0, B,[v, , v2.) = El Oca )) = M,, 
ax, 
(26) 


d(x 'y (L)) 


B,[v, , v2] = v,(L) = 0, B,[v, , v2] = “de - = @, 


It can be readily verified that the following expressions for v, and v, satisfy (25) and (26): 














BE sig, mw = Sy (ay 2 
Mi" a” ZL L°A \e, 3’ 
(27) 
EI Zi : x 
vo = — — : + = ’ 
ML 6L 2L 3 
where 
9 E : 
°o=- and G@ =k  . 
p p 
The normal modes uj and u3 which satisfy the equations’ 
i a (23'u" ae ; du’ me " 
EI — 2 ui) — uk As (az'u") + —*> = —plwi(xz'u5), 
dx; ( dx) 
(28) 
d(x5 M0") dur | @ 
rad an 3d 3? = —pAaw,us 
n | - df, + dx?! ‘iii 





‘In the remainder of this section the summation convention is suspended for the index n only. 
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and the boundary conditions 


u3(0) = uX(L) = 0 





(29) 
d(xz'ux(O)) _ d(xz'ur(L)) _ 0 
dx, ad dz, 7 
are 
xt2'u; = C, cos — ; u; = D, sin a . (30) 
Substitution of (30) into (28) yields 
c,| r(%") + pkA — plat | +D uta nr | 0, 
(31) 
2 
of | 4 p| (2) in (=) | - 0, 
L Co 
from which 
wy — | (wz) (i + 3) + Aah + (1") cics = 0, (32a) 


F A 2 4 1/2 
wo, = | (22) (ci + 2) + 4] oa {[ (x2) (ci +c) + 4a] - 4(%) ct} . (32b) 
There are two frequencies for each mode shape which correspond to the presence of two 
non-zero components of the displacement; they will be denoted by w,; and wp . 
If we return to (28) and (29), we may observe that these equations are also satisfied by 
z;'u, = C,, uy =0 (33) 


provided w, = A/Ic? ; u; is the so-called ‘‘thickness shear mode”. The notation is not 
meant to imply that the lowest frequency is wo ; in fact, this frequency is quite large. 
When (31) are solved for the ratio of the constants, the normal modes take their 


final form as 


aii nTx, ‘ (nwL)c; . NTL, 
2,'u; = cos —— “se = ss sin ——* 34 
vi . L , wl? — (nmc.)” L ( ) 


and the shear mode becomes x;'u} = COS wot. The solution will be complete when the 
integrals appearing in (14) and (17) are evaluated. Since f(é) is the unit step function 
a°f(t)/at? = 0, af(0)/at = 0, and f(0) = 1. Hence, P,(¢) = 0, and if we choose the initial 
conditions to be zero, the only non-vanishing constants are B, given by 


B, = -[ FE uw, + ux. | ae, / | E (ut)? + (ud) ‘| dz, . (35) 


Upon evaluating the above, we find after some manipulation 
2 2\-1] [@n2 j C; , 
2m — Ona) | (oo) — |r) | 
” my c,\’ 
ates, — o&){ (*2)' - (@)']. 


Buy 


(36) 


Bio 
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Finally, combining all of the above results, we have 


EI ry ra I (2) l Ie, 
y i << — —_-— —-——— COB Wot 





| i i te ie ee 3° L?Ac 
—2 i (war — wre) al A il (2) (<:)'] COS w,,f 
“ | (2x) = (2) | cos wath (37) 


This example has been solved previously by Leonard and Budiansky [8], using the 
method of Laplace transform. Equations (37) may be regarded as the “‘indicial ad- 
mittance” for this boundary driven problem. We could combine these results with 
Duhamel’s integral to obtain the solution for a much more general forcing function. 

5. Concluding remark. In the case where the non-homogeneous boundary con- 
ditions involve surface tractions only, it is possible to obtain a formal solution by using 
Lagrange’s equations of motion, such a solution is assumed to have the form given by 
(1) while the surface tractions are incorporated into q,(¢) through the generalized work. 
Thus the solution obtained in this manner will never satisfy the required boundary 
conditions of the problem, as may be seen from the fact that u? satisfy equations of the 
type specified by (8). 

Since the normal modes satisfy the homogeneous boundary conditions, it becomes 
necessary to establish the validity of the solution by showing, through a limiting process, 
that the solution satisfies the non-homogeneous boundary conditions as the boundary 
is approached from the interior of R. This process, quite often, requires that the solution 
also be obtained by another method, such as the one given in this paper. 

It may be mentioned that a similar situation may also arise in connection with the 
use of the integral transform methods when the solution is expressed in terms of a con- 


volution integral. 


REFERENCES 


. R. Courant and D. Hilbert, Methods of mathematical physics, Vol. I, Interscience, 1953 

2. R. D. Mindlin and L. E. Goodman, Beam vibrations with time-dependent boundary conditions, J. Appl. 
Mech. 17, 377-380 (1950) 

3. G. Herrmann, Forced motion of elastic rods, J. Appl. Mech. 21, 221-224, (1954) 

A. Clebsch, Théorie de l’élasticité des corps solides, translated from the German by B. de Saint-Venant 

and Flamant, Paris 1883, p. 129. Also see A. E. H. Love, Mathematical theory of elasticity, Dover, 

1944, p. 180 

5. R. D. Mindlin, Influence of rotary inertia and shear on flexural motions of isotropic elastic plates, J. 
Appl. Mech. 18, 31-38 (1951) 

6. R. D. Mindlin and H. Deresiewicz, Timoshenko’s shear coefficient for flexural vibrations of beams, 
Proc. Second U.S. Natl. Congr. Appl. Mech. 175-178 (1955) 

7. 8. P. Timoshenko, On the correction for shear in the differential equations for transverse vibrations of 
prismatic bars, Phil. Mag. ser. 6, 41, 744-746 (1921) 

8. R. W. Leonard and B. Budianski, On traveling waves in beams, NACA TN 2874 (1953) 


— 























51 


ON A GENERALIZATION OF TAYLOR’S VIRTUAL MASS 
RELATION FOR RANKINE BODIES* 


BY 
L. LANDWEBER 
Iowa Institute of Hydraulic Research, State University of Iowa 


1. Introduction. In 1928 G. I. Taylor [1] published a relation which expressed 
certain added mass coefficients of a moving body in terms of the singularity distributions 
which may be considered to generate the flow field about the body. This relation had 
the great advantage that, for Rankine bodies (i.e. for bodies generated by prescribed 
singularity distributions) these added masses could be computed directly from this 
relation without the necessity for the evaluation of the usually cumbersome expressions 
for the added masses as surface integrals in terms of potential functions. A derivation 
of Taylor’s formula has also been given by Lamb [2]. 

In 1950, in his Hydrodynamics [3], Garrett Birkhoff casually included a generalization 
of Taylor’s result. Recently, in the course of a seminar on added masses at the Taylor 
Model Basin, a new derivation of Birkhoff’s generalization was obtained by the present 
writer which, furthermore, resulted in a simple and elegant interpretation of a term 
which Birkhoff had left as a surface integral. These results will now be presented. 

2. Kinetic energy of a rigid body. Consider a rectangular Cartesian coordinate 
system (x, , 2 , £3) fixed to the body with origin at a point 0. The motion of the body can 
be described in terms of the vector velocity of translation U of 0 and rotation of the 
body with vector angular velocity w about 0. Let us denote the components of U by 
U; , U2 , Uz and the components of w by uw, , Us , Ue . Let o denote the mass density of 
the body and r the position vector of a point of the body with respect to 0. 

In vector notation the kinetic energy 7, of the body may be written in the form 


27, = / o(U +w X12) dr, (1) 


where dr is an element of volume of the body and the integration is taken over the 
volume of the body. By expressing (1) in terms of the components of the vectors, the 
kinetic energy can readily be expressed as a quadratic form in the velocity components [4] 
6 6 

oT", = y Muu; ; 

i=-1 -1 


P| 
where 


M;; ad M;; ’ 


M,, = M2. = M;; = M, the mass of the body, (2) 
M,; = My = My = My = Mos = Mae = 0,7 
Mi, = —M;s5 = Mri; Mu = —Mie = M23 ; 

M,, = —My = Mz, 


*Received March 28, 1955. 
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where x{ , x3 , x3 are the components of the center of mass of the body, My, , Mss, Moe 
are the moments of inertia of the body about the z, , x. , 2; axes respectively, M;, , 
Mes , M,; are the negatives of the cross products of inertia about the 2, , x. , 2; axes 


respectively, viz. 
M,, = -| o22x3 dr, ete. 


3. Kinetic energy of the fluid. It will be supposed that the body is moving through 
an inviscid, incompressible fluid, extending to infinity in all directions and at rest at 
infinity, and that the motion in the fluid is entirely due to that of the body. Then there 
exists a single-valued velocity potential ¢ which satisfies Laplace’s equation 


V%=0 (3) 
and the boundary condition on the surface of the body 
dp _ 

dn CrwxXH >a, (4) 


where n is the unit vector normal to the surface of the body, directed into the fluid. Let 
N; , M2 , M3 denote the direction cosines of the normal. Then the boundary condition (4) 


may also be written in the form 


6 
do is 
= n;U d) 
dn > sate ( 
where 
Ne = XWoN3g — XM3Ne ’ Ns = M3, — Nz ’ Ne = GN2 — Gl ’ (6) 


i.€. M4 , Ns , N. are the components of the vector r X n. 
Because of the linearity of Laplace’s equation and the boundary conditions, we 
may now make the Kirchhoff assumption that 


6 


$= dud. (7) 


t=1 


Here ¢, is the velocity potential when the body has only the one component of motion, 
u; , and that of unit magnitude. From (5) and (7) we obtain the boundary conditions 


de, 


“atta (8) 


Let 7, denote the kinetic energy of the fluid and p the density of the fluid, assumed 
to be uniform. Then we have 


ee db 49 
ar, = —p |[ ¢ as, (9) 
where the integration is taken over the surface of the body. From (5) and (7) then 


6 
2T. = i: Aj ju; ’ (10) 


6 
i=l 
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where 
Au =p |f 6, a8 (11) 
are the added masses of the fluid. It may also be shown that 
Ai; = Aj; . (12) 


The foregoing brief review of the theory of added masses is essentially a summary 
of the treatment in Lamb’s Hydrodynamics [5). 

4. Taylor’s formula and Birkhoff’s generalization. ‘Taylor’s formula may be written 
in the form 


Aj; + M - dor py; ’ a = 1, 2; or 3, 


where u;, is the 7th component of the total dipole moment of the sources and doublets 
which give the flow about a body translating in the 7th direction, and M is the mass of 
the displaced fluid. Birkhoff’s generalization [3], in our present nomenclature, becomes, 


Ai — P II x; ie dS = 4rpyi; , t= 1,2,or3;j = 1,2, --- , 6, 
where u,; is the 7th component of the dipole moment of a body corresponding to a j-com- 
ponent of motion of translation or rotation, and the surface integral is taken over the 
surface of the body. 

Birkhoff derives his generalization by applying Green’s reciprocal theorem to 
x,(do;/dn) — ,;(dx;/dn) for a region bounded internally by the body and externally 
by the surface of a large sphere. Lamb [5] derives Taylor’s formula by equating 
asymptotic expressions for the potential, given first in terms of the added mass, and 
then in terms of the dipole moment. It is shown in the following sections that Lamb’s 
method can also yield Birkhoff’s generalization, and furthermore that the term con- 
taining the surface integral is simply M,, , the corresponding component of the mass 
tensor (2), foro = p. 

5. Asymptotic value of potential at a great distance from body. Let r denote the 
position vector of a point on the surface of the body with coordinates x, , x, , Z3 , and 
R the position vector of a point P at a great distance from the body, with coordinates 
£, , & , &; . Let R be the magnitude of R and R, the magnitude of the vector R — r, i.e. 
the distance between P and a point on the surface. 

The potential at P is given in terms of the values of the potential and its normal 
derivative on the surface by the well-known formula [6] 


r a(i -}4| ; 
—~ // E dn (z) R, dn ” ” 


Let us expand 1/R, in a Taylor series. We obtain 


3 ee & Dib + Xb. + tabs ah (14) 


R, R R® 





and, from (14) 


d {1 0 0 re] 
dn (4) ™ (n on, * ae, * ™ Se 





1 oo Miki + Nak. + Nags pues 
Fs R " - 
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Then, from (7), (11), and (15) 


ft, = 
J] 6% (z)4 dS ~ Fp i] yu % 2 > ne, ds 


or 


1 6 3 . . 
[fea ? ain L(t) as ~ OR & LA Ub; - (16) 


Also, from (14) 


ff -1#as~ ff fe it tab + ts) d¢ ag 
J. R, dn ° t dn 


) since the total flux of velocity through a closed stream surface 








But {ff — (d¢/dn)dS = 
is zero. Hence, from nt 


I - - de 5 a = | Ye wn | Y 2g, dS ~ ds DD 3 or uk: [[ ne, dS. (17) 





ity an” i 
But for 7 = 1, 2, or 3, Gauss’ transformation gives 
ss Ox; ene 
I n,;dS = I i dr = 6,;V, (18) 


where 6,; is the Kronecker delta and V is the volume of the body. When z = 4, however, 


we obtain from (6) 


| nt; 48 = I (x2%j;N3 — %3XjN2) dS, 


ve 


[ if j=1 
= [ff (eed. — 295,:) de =} —23V j=2, (19) 
ey fas 


where x{ , x3 , x; are the components of the position vector of the centroid of the volume 
of the body. Comparison of (18) and (19) with the values of M,; tabulated following 
Eq. (2) shows that, when the body is of neutral density (¢ = p), 


|| nia. c1%- 4360648 (20) 
Hence, substituting the results in (16), (17), and (20) into (13), we obtain 


4rpR*op = > y (Ai; + M,,)ud; . (21) 


t=1 


Now suppose that the flow due to unit value of the 7th component of the motion is 
generated by a set of sources C,,; and a set of vector doublets d,; with components 
das » Ania » Anis » Then the potential at a point P due to this system of sources and doublets 


is 
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6 


= Eu (Det Dm) cm 


t=1 


where R,,; is the distance from the mth source to P, R,,; the position vector of P relative 
to the nth doublet and R,,; the magnitude of R,; . But, from (14) 





ates is L + Lmikti + Lmok2 + mats 
a R® ’ 
where 2; , mo» Um3 are the coordinates of the mth source. Also 
Cai po 1 = 
UR RY Cm = 0 


since the totality of sources generating a closed body is zero. Furthermore the last term 


in (22) becomes asymptotically 


—— 





oe oe 2 
} od ~~ — > P 
Hence we obtain from (22) 
6 3 
R*¢p = p i (> Cuitai + } d,,;;)Uié; . (23) 
i=l j=l m n 


6. Relations between added masses and singularity distributions. Finally, com- 
paring the asymptotic values of the potential in (21) and (23) we obtain the desired 
relation between the added masses and the source and doublet distribution 


Ay; + Mi; = 440( >> Caitni + D> dais) (24) 


¢= 1.2. --- 63 = 1, 2,38. 

It should be noted that Eq. (24) does not give relations for all the added masses. 
The six terms due to pure rotation, such as Ay, , Ay; , are not given. Attempts to find 
such relations, either by taking additional terms in the asymptotic series for the potential 
or by using Birkhoff’s method [3], have given numerous results which resemble but are 
not quite the ones sought. Thus, although a relation for Ay, = ff @s(r2n3 — 23n.)dS 
has not been found, it can be shown that 


I 4(LoNzg + TgN2) dS + If (x3 — x3) dr 


== 4n[ >> CoMasbus + > ie (daso%ns + dna3Tn2)], (25) 


where the volume integral is taken over the volume of the body and z,, , Zn2, Zn3 are the 
coordinates of the nth doublet. It is interesting to compare the term My = JJf 
(x3 + 2x3)dr, which would occur in Eq. (24) if it could be applied for 7 = 4, with the 
corresponding term in (25). Also it has been found that, although a relation for Ay; = Jf 
$;(X2N3 — 23N2)dS has not been found, we have 


// $s(r2Ns + xn2) dS — Ii 2%, ar 


= 4n[ >> Cuslm2tm3 + > (dys2tN3 + dassTn2)]. (26) 
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In this case the term M,,; = — {ff x,2.dr occurs explicitly. These relations suggest certain 
guesses as to how (24) may be generalized, but we will refrain from mentioning them 
until some evidence for their validity is adduced. 
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A PRIORI BOUNDS FOR TEMPERATURE IN 
CIRCULATING FUEL REACTORS* 


BY 
JOEL FRANKLIN anp HERBERT B. KELLER 
Institute of Mathematical Sciences, New York University 


1. Introduction. The neutron density, u, and temperature, v, in a circulating fuel 
reactor satisfy non-linear partial differential equations of the form 


U, = Uzz + Alvju, (1.1) 
v, +o, =u O<zx<a,t>0). 


These equations are in dimensionless form; ¢ is a time-variable; z is a space-variable; 
and ¢ is the positive, constant speed at which the fuel flows through the reactor. A (v) 
is a given function of v. The neutron density is zero at the boundaries of the reactor: 


u(0, t) = u(a, t) = 0. (1.2) 
The input fuel is kept at a constant temperature, say 
v(0, t) = 0, (1.3) 
in an appropriately chosen scale. The initial state of the reactor is given by 
u(x, 0) = u(z), v(x, 0) = v(x), (1.4) 


where w% , Y¥ are known functions. A physical derivation of the above equations, with 
A(v) = constant-v, appears in a forthcoming report by J. Fleck [5]. 

The purpose of this paper is to determine a priori bounds for norms of wu and v. 
Under natural assumptions on the given functions A, Uo , vo , we suppose that smooth 
solutions u, v exist to the mixed initial- and boundary-value problem posed by the 
preceding equations. We then derive numerical bounds, depending only on the given 
functions, for the maximum-norm of v and for an integral-norm of u. These bounds will 
hold in the infinite domain 0 < x < a,0 < t¢ < o. Our results are found by the elemen- 
tary sort of argument used by E. Hopf [1] and by L. Nirenberg [2] in their justifications 
for the maximum principles for elliptic and for parabolic equations. 

To prove a preliminary result, we shall use the following form of the weak maximum 
principle for parabolic equations. Let U(z, t) be a solution of the equation 


U, = U,, + ¢(2, DU (1.5) 
which is twice continuously differentiable in the closed rectangle 


R;: 0<2zs4, 0O<t<T; (1.6) 
*Received April 1, 1955. This work was sponsored by the United States Atomic Energy Commission. 
The authors wish to acknowledge their indebtedness to J. Fleck of the Brookhaven National Laboratory, 
who suggested a specific problem from which the present work originates. 
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let o(x, #) be continuous and < 0 in R;, ; assume that U takes positive values in R, ; 
then U attains its maximum value on one of the boundary segments 


x= 0, { = 0, 2 =a; (1.7) 


similarly, if U takes negative values in R; , then U attains its minimum value on (1.7). 
This principle has long been known in various forms; see for example M. Picone [3]. 
In [2] Nirenberg proves the strong maximum principle, which asserts that, if U attains 
a positive maximum or a negative minimum at any point of R, not in the set (1.7), then 
U = constant. 

2. Assumptions on the given functions. We suppose that the given functions 
Uo(X), Yo(x) are twice continuously differentiable* functions of x satisfying the boundary 


conditions 
u(0) = u(a) = 0, v,(0) = 0. (2.1) 
Because u represents a density, we shall require 
u(x) > 0 O< 2 <a). Gz 
Let 
a = min»,(2), B = max v,(z) O<2z< a). (2.3) 


We suppose that A(w) is a given, continuous function defined for all w > a. Since a 


< v(0) = 0, we may define 
Bw) = | Aw) dw — (w > a). (2.4) 
We now make the important assumption 
lim sup B(w) < min (7, 0), (2.5) 
where 
y = min [—w(x) + cvi(x) + v6’(x) + BO,(2))] O<2<a). (2.6) 
This assumption does not appear to be unnatural; for example, if A(w) = constant-w, 


we have B(w) — — © as w —~-, so that (2.5) holds regardless of the value of y. Since 
B(w) is a continuous function, with B(O) = 0, assumption (2.5) implies that there exists 
a unique, finite, non-negative number M such that 

9 


B(M) = min (7, 0), and B(w) < min (7, 0) for w> M. (2.7) 


3. Boundedness of temperature. 
Theorem. Let u(z, t), v(x, t) be defined for 0 < « < a,0 < t < © as functions with 


continuous derivatives up to the third order. Let u, v satisfy the differential equations 
(1.1) and the boundary and initial conditions (1.2), (1.3), (1.4). Let the given functions 
Up , Vo , A satisfy the conditions of Sec. 2. Then v(z, t) is bounded; v satisfies the inequality 

a < v(x, t) < max (8, M), (3.1) 


where a, 6, M are the constants defined in (2.3) and (2.7). 
*We will not attempt here or in the statement of the theorem to use the weakest possible assump- 


tions of smoothness. 
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Proof. As a preliminary result, we prove that u(z, t) is > 0 for all ¢ > 0. Consider 
the rectangle R, defined by (1.6). Following Nirenberg, we define 


\ = X(T) = max Alp(z, 0d] for (xz, i) in Rr 
and define 
U(x, t) = eu(z, 2). (3.2) 
Then U satisfies Eq. (1.5) with 
o(x, t) = Afo(z, )] — d. 


Since ¢ is < 0 in R; , we may conclude from the maximum principle that U, and there- 
fore u, assume negative values in R, only if U assumes a negative minimum on one of 
the segments (1.7). But u, and therefore U, are > 0 on (1.7) because of the boundary 
conditions (1.2) and the inequality (2.2). Therefore, u is > 0 in R; . Since T is arbitrary, 
u(x, t) is > Oin R, , i.e. for all t > 0. 

To prove the theorem it suffices to show that v(z, 2) satisfies the inequalities (3.1) 
in every finite, closed rectangle R; . We consider RP; as the sum of four boundary segments 
and the interior, according to the following definitions: 


N: 0<2z <4, t=T; (3.3) 
E: Z= 4, 0<t<T; (3.4) 
S: 0<-2 <a, t= 0; (3.5) 
W: z= 0, S<tg?; (3.6) 
i: 0<2<e, e<i<T. (3.7) 


We also define sets R7 and R7 as the intersections of R; with the half-planes ct — z < 0 
and ct — x > 0. 

To find a lower bound for v, we draw from each point p in R; the segment of positive 
slope 1/c extending to a point p’ on S + W. The point p’ lies on S or on W according 
as the point p lies in R; or in R; . Since u > 0, it follows by integration of the equation 
v, + cv, = u that v(p’) < v(p). By (2.3) and (1.3) we have v > a on S and v = 0 on 
W. Therefore, 


v>ainR; and v>OimR;, (3.8) 
from which, since a < v(0) = 0, we find the required lower bound 
v > a in Rr . (3.9) 


Using the definition (2.4) of B, we may eliminate u from (1.1) to obtain the single 
equation for v 


0 0 
This equation implies that there exists a continuous function f for which 


v, = v,, + BY) + fect — 2). (3.11) 
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To find an upper bound for v we shall need an upper bound for f. On S (3.11) becomes 


f(—z) = o,(z, 0) — v,,(2, 0) — Biv(z, 0)], 
= U(x) — cvi(x) — v6’(x) — Blvo(x)]. 
Therefore, by the definition (2.6), 


fy <-y (-a<y<0O). (3.12) 
On W (3.11) becomes 
f(ct) = v,(0, ) - 


v,(0, t) — Biv(0, a]. 
But v = v, = 0 on W, and B(0) 


v, = 0 on W. But, by (3.8), v(p) 


(3.13) 
> 
point of W. Therefore, 


0. Since u = 0 on W, the second equation (1.1) gives 
Oin R;, 


so that v(p) > 0 in a neighborhood of each 


v., > Oon W. 
It now follows from (3.13) that 


(3.14) 
fly) < 0 (y > 0). (3.15) 
Combining (3.12) and (3.15), we find 
. ; —¥ for 
flict —xz< 


x,t) mR; ' ' 
(, : (3.16) 
0 for (x,t) nR;. 

Let P be a point at which v attains its maximum value in R; . Then P lies in exactly 
one of the four continua W, S, N + I, E. If P e W, then v(P) = 0, by the boundary 
condition (1.3). If P e S, then v(P) = 8, by (2.3). 


Suppose P e N + J. Since v(P) is the maximum value of v, and since the set N + I 
is open with respect to z and open from below with respect to ¢, we must have 


oF) = 0, 


v.(P) < 0, 
Therefore, by (3.11) and (3.16), 


v(P) > 0. (3.17) 


Bi(P)] > —f => min (7, 0). 
From (2.7) we now conclude 


(3.18) 
uP) < M. 


(3.19) 
Finally, suppose P e E. Since v(P) is the maximum value of v, the definition (3.4) 
of E leads to the inequalities 


v(P) > 0, v(P) > 0. 


(3.20) 
But v, + cv, = u = Oon E, and cis > 0. Therefore, (3.20) yields 


v(P) = 0, v(P) = 0. 
From the maximum-property of v(P) we may now conclude 


v.AP) < 0. 
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From (3.21) and (3.22) follow (3.17) and, hence, the inequality o(P) < M. Summarizing 
the results for all of the four sets W, S, N + J, and E, we find v(P) < max (8, M) re- 
gardless of where the maximum value v(P) is attained in R; . This completes the proof 
of the theorem. 

4. Boundedness of a norm of neutron density. By comparing the equations 

u, = Uz, + Alvju, ut = ux + ru*, 
where \ is a constant > A(v), one easily obtains (see Pélya and Szegé [4]) an inequality 
of the form 
u(x, t) < K exp (A — x°/a’)t, (4.1) 
where K is a positive constant. According to (3.1) we shall certainly have \ > A(v) 
if we define 
A = max A(w) la < w < max (8, M)]. 

Then (4.1) shows that the non-negative function u is bounded if the length a is suf- 
ficiently small, i.e. if a” < ’/\. But this restriction is too strong in cases of physical 
interest. 

If a is unrestricted, the arguments of Sec. 3 evidently do not yield an upper bound 


for the neutron density u(x, t). However, as an immediate consequence of the theorem, 
we can find a bound for a certain integral-norm of u. Define the boundaries 


E.: z=a, t>O0; S: O0<2z<a, t=0; W.: x=0, t¢>0. 


] 


From a point on Z., , with ordinate ¢, we draw the segment L, of slope 1/c to a point on 
S + W.. We then define 


n() = [ u ds; (4.2) 


this integral may be properly called a norm, since u is > 0. By the second equation 
(1.1) we have 


dv 
—=uonL, . 
ds 


(1 + > ae 
Therefore, 


n(t) = (1 + c’)'”’ (a, ) — v(a’, ¢’)), 


where (a’, t’) equals (a — ct, 0) or (0, ¢ — a/c) according as ct — a < Oorct —a>0O. 
Thus, by (2.3) and (1.3), 


> anlis a if ct-—a<0 
—v(a’, t’) 
= 0 if ct—a> 0. 


Since, by the theorem, v(a, t) < max (8, M), we have proved the following result: 
Corollary. Let the conditions of the theorem be satisfied, and let the norm n(¢) be 

defined by (4.2). Then 

max (8, M) — a if ct—a<0 


max (8, M) if ct—a>O0O0. 


(l+e¢)"'’n(t) < 
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THE DIFFRACTION OF A CYLINDRICAL PULSE BY A HALF-PLANE* 


BY 
ROBERT D. TURNER (Ithaca, New York) 


I. Introduction. In the past fifty-eight years, many papers have been written on 
the theory of diffraction in two dimensions by a wedge. Since the time of Sommerfeld’s 
attack on the half-plane [1]**, several approaches, using various specialized forms of 
excitation, have been investigated. Sommerfeld’s use of multi-valued functions for 
plane wave excitation has been supplemented by Macdonald’s expansion in the appropri- 
ate eigenfunctions for the problem [2]. Much of the work that has been done in this field 
has been for the domain of harmonic time-dependence; in a sense, these results represent 
complete solutions for the particular form of excitation involved. However, the inverses 
of these Fourier-transform solutions are rarely (if ever) determined. 

Recently, more direct attacks on the time-dependent wave equations have been 
successful. Keller and Blank [3] have obtained useful results for the scattering of plane 
waves by wedges and corners. By recognizing the nature of propagation of discontinuities 
for solutions of hyperbolic partial differential equations, they are able to make an 
appropriate change of independent variables. A change of dependent variable yields 
the Laplace equation, which is then solved by conformal mapping. More recently, Kay 
[4] has achieved a rather general result, in that he is able, at least in principle, to write 
down the solution for an arbitrary form of excitation. The method is to make a change 
of independent variables suggested by the work of Keller and Blank. This leads to a 
non-orthogonal co-ordinate system, and he devises an integral transformation to treat 
the rather complicated differential equation he obtains upon separation of variables. 
The transform kernel involves a Whittaker function, and the integration is over a 
reasonably complicated contour in the complex plane of the separation parameter. 
This brings up what seems to be the only limitation on his method—the mathematician’s 
ingenuity in carrying out the details of the analysis. The intricacy of the transformation 
precludes the possibility of obtaining a more general result than he does—viz., the 
verification of the results of Keller and Blank. 

The method of attack that we propose is less sophisticated than the foregoing, but 
possesses the advantage that we can determine the Green’s function for the problem. 
For simplicity, we restrict ourselves to the case where the wave-function is constrained 
to vanish on a half-plane, the excitation being a line source parallel to the edge of the 
half-plane. The analysis is no more difficult for the general case, however, and may 
even be extended to the diffusion equation. 

II. Explicit statement of the problem. We wish to determine a scalar function 
@ such that 

1d dr — r) 5(8 — H) H(t — 0°) 


al las c at To (1) 





subject to the boundary condition 
¢=0 (2) 
*Received April 6, 1955. This paper is part of a thesis submitted in partial fulfillment of the require- 


ments for the degree of Doctor of Philosophy in Applied Mathematics at Harvard University. 
**Numbers in square brackets refer to the bibliography at the end of this paper. 
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for 6 = O and 6 = 27 and the initial conditions 


- ; Op _ 
@¢=0 and Po = 0 (3) 


at 
t = 0. 


Figure 1 shows the relationship between the source point (ro , 0), the observation point 







(8) R 
OBSERVER (6, 0.) 


SOURCE 


OBSERVER ON SCREEN 








SCREEN 
Fig. 1. Relationship between source, observer, and diffracting screen. 
(r, 6), and the diffracting screen 6 = 0 (or 6 = 2"). For convenience, we denote by R 
the distance from the source point to the observation point: 
R = [r’? +175 — 2rr, cos (@ — @)]'”’; (4) 


when the observation point is on the screen, we call the distance R, : 


R, = [z? + ri — 2zro cos ]'”’. (5) 
The problem is reduced to a simpler boundary-value problem at once by writing 
¢=¢o"+4u (6) 
where [5, p. 332, Eq. 57] 
inc 0 t < R/c ~ 
¢ a (7) 


((1/2m[@ —R/e}? t>R/e 
¢'”° is the free-space Green’s function for the two-dimensional wave equation, and 
satisfies the inhomogeneous wave equation (1) as well as the initial conditions (3). 


Accordingly, 


' 1 3’ 
Vu = : = (8) 


and on the screen (@ = 0 and 6 = 2n), 
/ 
ee t < Ro/e 0) 
((—1/2m)(#@ — R2/e)"'? -t > Rife 


Also, u satisfies the initial conditions (3). The problem is then to find wu. 
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III. Solution of the problem. We now consider the one-sided Laplace transform 
of u, which we denote by U. On making use of the initial conditions (3), the partial 
differential equation transforms to 

V°U — (p’/e)U = 0, 
and the boundary condition becomes [6, p. 125] 


U = (—1/2n)Ko(pRo/o) 


on the screen. The function K, is the zero-order Macdonald function [7, p. 78], which 
is proportional to the zero-order Hankel function of the first kind with imaginary 


argument. 
It is convenient to discuss the function 


v= U + (1 /2n)K(pro/e), (10) 
which vanishes with r. After multiplying by r’, the partial differential equation satisfied 
by v is: 

roy roy ov pr _ PT oy / 
or” + or + 00° line c (1/2m)Ko(pro/c) (11) 


and the boundary condition on v is that 
v = (1/2n)[K.(pr./c) — Ko(pR/c)] at ¢6=0 and 6 = 2x. (12) 
We now make use of the Griinberg modification [8] of the integral transformation 
of Kontorowich and Lebedev [9, 10]. This states that if 
V(s) = | v(pr/c)K (pr/c) dr/r, (13) 
where K, is the Macdonald function of order s, then 


(pr/) = (1/ri) |  V@.(pr/e)s ds, (14) 


where J, is the modified Bessel function of order s [7, p. 77]. In transforming (11) by 
means of (13), we integrate by parts several times and make use of the modified Bessel 


equation satisfied by AK, , obtaining finally: 


OV, oy, _ _ 8Kolpro/c)_ . 
ag + ©! = —Tain (wa/2) (15) 


Here we have made use of the formula 


ro Ts 
[ Kis és = se 


[7, p. 388, Eq. 8; 6, p. 1]. The boundary condition (12) must also be transformed; at 
6 = O and 6 = 27, we have V = V,, where 


V, = (1/2) / K.(pr/c)(Ko(pro/e) — Ke(pRo/e)] ar/r. (16) 
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This integral is evaluated in the appendix. We have now reduced the partial differential 
equation to a particularly simple ordinary differential equation. All that remains is to 
solve the ordinary differential equation and carry out the two inversions for the two 
transformations used. The solution to the differential equation is 


4 / 
V(s) = A(s) sin 6s + B(s) cos 6s — - Ko(pro, 2... 
4s sin (1s/2) 





On applying the boundary condition, we have 


Kolo | cos (x — 6)s Ko(pro/c) ; 


4s-sin (28/2) COS 78 ~ 48-sin (xs/2) 


Vis) = lv. + 
Inserting the expression for V, , we have: 


cos (r — 4)8s cos (r — @)s K.(pro/c) 


8 sin 2s 4s sin (xs/2) 





V(s) = K,(pro/c) 


Note that V is analytic at s = 0. (This is easily seen if we recall that K,(z) is even in 
s [3, p. 79, Eq. 8], so that for s > 0, K,(z) = Ko(z) (1 + as’ + ---).) 
The next task is to invert this transform: 


tf | K.lor g Se Shoe oS [tories ds (17) 


“= ; : 
Tt Jie 8 sin 27s 4s sin (xs/2) 





For r < 7) , we may use the inversion formula (14) as it stands, computing the integral 
by closing the contour with a semicircular arc on the right as in Fig. 2, and applying 


5 - PLANE 











Fic. 2. Contour for inversion of the transform. 


the theory of residues. If we choose the radius of the semicircle to be } (2m + 1), where 
m is @ positive integer, it is not too hard to show that the integral over the semi-circular 


arc must vanish as m ©. Thus, 


vy = —2 )- residues inside the contour = v, +», , 


where 
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o, = —2 >" residues of K,(pr./c)I(pr/ec) COs (x — 6)s cos (xr — 6)8 
nol sin 27s 
ats = n/2; 
— K.(pro/c)I.(pr/c 
v, = —2 > residues of hae fon ) t s = 2n. 


We evaluate »v, first, since it is the simplest: 
v2 = (1/m)Ko(pro/e) 2) (—1)"Taa(pr/c). 
n=1 
This is easily summed by means of a well-known summation formula for Bessel functions 
(7, p. 34, p. 77]: 
v, = (1/2n)K.(pro/c) — (1/2m)Ko(pro/c)Io(pr/c). 
Note that the first term of v, is exactly what we added to U in Eq. (10). Next, 


v, = (-—1/n) =. (—1)"Ky/2(pro/c)I,/2(pr/c) cos 5 n( — 6) cos + n(x — 6). 


n=l 


After some trigonometric manipulation, we find that 


v = (1/2n)Ko(pro/e) + (1/m) D2 Kuya(pro/e)Ina(pr/c) sin t né@ sin , nO, 
n=1,odd 2 2 
‘| = | 
~ {Kotor c)I,(pr/e) + 2 >> K,(pro/o)I,(pr/c) cos n(@ + ay} 
n=1 


ae {Kater c)I(pr/e) + 2 >> K,(pro/c)I,(pr/c) cos n(@ — aa}. 


dr n=1 
The bracketed terms may be summed [7, p. 361, Eq. 8], giving finally: 
v = (1/2x)K,(pro/c) + (—1/42)[Ko(pR/c) + Ko(pR,/o)], 


+ (1/m) SS Karlpro/c)In/2(pr/c) sin ; né@ sin 5 nO , 
dd 


n=1l,o 
where 


R, = [r? +15 — 2rro cos (6 + 6)]'”; (18) 
R, is the distance of the observation point from the image (with respect to the plane 
y = 0) of the source point. Finally, 
: , — 1, : 1 ’ 1 
U = (1/r) >> Kasisolpro/0)Ins12Apr/c) sin (n + 5/9 sin \n + 5 A 


— (1/4)[Ko(pR/c) + Ko(pR,/c)]. (19) 


Equation (19) is valid only if r < ro. If r > ro, it is necessary to revise the form of 
the inversion integral. This is done by utilizing the defining equation for the Macdonald 
function [7, p. 78, Eq. 6] and the symmetry of V(s), which is always even. Doing this, 
Eq. (14) becomes 


v = (1/m1) [ V(s)(—2/x)K,(pr/c)s sin xs ds. 
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> 
@ 


Now when we insert the expression for V(s) into the inversion formula, the term in- 
volving K,(pr,/c) will give the same result as before, when we sum over its residues. 
The other term turns out to be just the integral for V, that we evaluated previously, 
with r and r, interchanged. This merely interchanges r and rp in the expression (19) 
for U. Introducing the shorthand notation: r, = max (r, ro), r- = min (7, To), the two 
expressions for U may be combined: 


U = (1/n) z. K,41/2(prs/e)In+1/2(pr-/c) sin (n + Ne sin (n + V6, 
— (1/4m)[Ko(pR/c) + Ko(pR,/c)]. (20) 


The next step is the inversion of U. This is facilitated by use of the integral of Sonine 
and Gegenbauer [7, p. 367, Eq. 17]. Specialized to our problem, this reads: 
Ko1/e(pr>/0)1n+1/2(pr-/o) 
exp |—(p/e)(r" + 79 — 2rro cos >) *~ ‘ 
= = (rr,)'” [ p | u SS +r ee $) sin ¢ d¢, 
, (r° + To — 2rro COS d) 


- J0 


when P,, denotes the Legendre polynomial of order n. But this last integral is the Laplace 


transform of 


Rio age [ 6[t — (1/c)(r* + 75 — 2rro cos ¢)'’*] ; 
(rr,)'/? ——, <3 P.,,(cos ¢) sin ¢ dd. 
(r? + r5 — 2rro cos ¢)'”” nCOS @) ita 


Performing the indicated integration, we find that 
Kas1/2(pr>/0)1q+1/2(pr</c) 


is the Laplace transform of the function 


~wypif trm—¢et 
~C(rTo) pl = ‘ | ir —% | < cl <P + fo 


(1 
42 217 
LO otherwise. 
The remaining terms of (20) are easily inverted [6, p. 125], so that u = uw, uz; — uw, 
where 
f ¢ r+rn—ctl | . ( La ( l 
| Qa(rro)”” > al Orr | centh ‘i + 3) 0 sin sited 2 4% 
- (21 
med r—nm|<e<r+n 
0 otherwise, 
0 t<R/c 
U2 = 4 | (22) 
| ——, -——{; t> R/e 
4n[t? — R°/c’] - ; 
[ 0 t< R,/e 
Us, = 4 a (23) 
t eo R, le. 
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Now ¢ = ¢'”° + u, and on comparison of (22) and (7), we see that ¢'"° = 2 uz , so that 
d= uy + U2 + Us. (24) 


This completes the formal solution of the problem. Since the series converges very 
slowly, it is not a very practical result as it stands. If the Green’s function @ were to be 
used for the determination of the response to some less singular excitation, the series 
resulting from term-by-term integration of (2) would converge much more rapidly. 
In special cases, however, the series (21) may be summed, and we shall now consider 
some examples where this may be accomplished. 

There is an apparent difficulty that may be resolved at once. The term u; represents 
a source at the image of the source point whereas there is no such source in the original 
problem. There is, however, a contribution from u, which just cancels the singularity 
in U, at (ro , 27 — 4). To see this, we investigate the behavior of u, and u, at this point. 


We have at once that 


uz, = — 1/4rt i> 0, 
and 
(c/2nro) > Pj 1— ef sin? {n + : 6 0 <ct < 2r, 
u, = J aTT > a n or b 5 fi} “lo 5 
0 otherwise. 
Using sin” x ‘(1 — cos 2x), u,; may be broken up into two series. The series involving 


cos (2n + 1)@ is not singular as t — 0. We therefore consider only 


if ' c Ce 
iit 4 Sorry) »» pi _ 3) 0 < ct < WM 
| 0 otherwise. 


This series may be summed [6, p. 53]: 
ul = 1/4at, 0<ct <2. 
Thus, uw; , as stated above, cancels the singularity of u; at the image of the source. 
Similarly, at (ro , 9), the u, term adds to the term u, to give the proper source-like 
character (i.e. with a coefficient 1/27 instead of 1/47) at this point. 
The next case to be considered is where 6 = 37/2, 0) = 2/2 (see Fig. 3). 





dy SOURCE 

| 

SCREEN 

; 

| 

| 

| 

(%G,37/2) ' IMAGE OF SOURCE 

| 


(% 37/2) | OBSERVER 


Fic. 3. Configuration for a test of Babinet’s principle. 
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The series may be summed as in the previous example, yielding 
g — \A/4M TE —@ +n) — [P-—G—n)/ey'") ct >rt+r 
{ 0 otherwise 
Now the problem complementary to the case under consideration has as the boundary 
condition on the complementary screen, y = 0, x < 0, 
0 
A = 0 
on 
The solution here is 


‘ > 


y= ((1/4m){((? — (r+ n) ey? + [@ — ¢ — 1)?/e)"'7} ct >r+r 
lo otherwise. 


Now Babinet’s principle states that the signal and the signal for the complementary 
problem must add up to the incident excitation. Adding our two solutions, we have 
(1/2m)[(? —(r+n)’/e)'? ct>r+n 
eres [ /c*] + 


lo otherwise, 


which is indeed the incident field. Thus, Babinet’s principle has been verified in the 
time domain for this special case. Because of the highly discontinuous nature of the 
solutions, this constitutes a rather severe test of the validity of Babinet’s principle in 
the time domain. 

Finally, we shall compute the “charge density’’—that is, the discontinuity in d¢/dn 
across the screen, for the case @ = 2. From (20), the Laplace transform of p(x), the 


charge density, is given by 


; =~ «. y F l 
P(x) = (2/x2z) > Kasso(prs/c)In+1/2(pr-/c)(— 1) (n + i). 
None of the other terms contribute to P(x) since their normal derivatives are continuous 
across the screen. This series may be summed [7, p. 366, Eq. 11], giving 
(1/n)(ro/x)'* exp [—(p/c)(x + 1)] 


P(x) = 
v) fot zt 


so that 
5[t —_ (1 c)(x + To)] 


we—— is .(r,/x)'”*. 
a) (ro + 2) »/ Z) 


” 


Thus, the charge on the screen has no after-effect, and the “amplitude” of the pulse 


falls off as 2~*”” as we move away from the edge. 

IV. Conclusion. We have derived Green’s function for the time-dependent wave 
equation and considered a few of its properties. The method is quite straight forward— 
two transformations and two inversions—and is readily generalized to the case of the 
wedge. As we remarked in the introduction, we may also treat the problem of (say) 
heat flow due to an instantaneous source near a wedge held at constant temperature or 
near an insulated wedge. The principal change to be made is the substitution of p'”’ for 
p in all of the equations up to Eq. (20). The inversion of the product 


| wap r le) Ins vip r/c) 
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is not as easy, but can be carried out for any specific n. Although the general term in the 
series is difficult to determine, the series converges rapidly enough so that this is not im- 
portant. (Note that the extremely slow convergence of the series (21) is not a fault of 
the method, but only a consequence of the character of the solution). A paper on this 


is planned in the near future. 
Appendix. Evaluation of V,. With a convenient change in notation, we have 


Vo = (1/2n) [ K,(x)[Ko(to) — Ko{(x? + 23 — 2x2» cos 0o)'7}] dx/zx. 
#0 
We use the addition theorem for Macdonald functions [7, p. 361, Eq. 8], and get 


2rV, = [ K(a)| Kaz — y €mE m(Xo)1_(Z) COs mo, | a 


m=O 


+ [ K.(a)| Kaz) nei > €mK m()I m (Xo) _ mo. 7 ' 


m=0 


where ¢,, = 1, m = 0;m = 2,m > 1. Interchanging summation and integration we have: 


az 


2V. = K,(2) | , K,(x)[{1 — I,(x)] = 2 >> K,,(to) cos m4 [- K,(x)I,,.(x) ar 


zr 


+ [ K,(x)[Ko(%o) — Ko(x)Io(x0)] a — 2 > 1,(20) cos méy | K,(x)K,,(x) af. 


vz 


The second and fourth integrals are evaluated by means of the indefinite integral 
for cylinder functions [7, p. 134, Eq. 7]; when this is done, the two series may be combined 


to give 


~27, 5 8% KK! — K1-K! — 1K.Ké + 1.KK'l, 


2 2 
noi 7 — & 





where the argument 2, has been suppressed. The second and fourth terms in the square 
brackets cancel, leaving just K,(z») times the Wronskian of K,, and J,, , which is just 
1/2, [7, p. 80, Eq. 19]. Thus, the series in the expression for V, combine to give 


=> £08_mB, 


—2K.(%.) 20 aa 


m=1 s 


This last series may be summed [11, p. 278, Eq. 13] so that the expression for V) now 
reads 


QV, = KolXo) [ K,(x){1 — In(x)] e — Ko(2o) [ K,(zx)[1 — I(x] “ 


-7\| CoS (er — O)s 1 [ dx ms dz 
+ rK,(29)| — 4] + Ko(xo) a K,(x) - To(2o) [ K,()Ko(z) > 
cos (x — %)s _ 4, 

8 sin 7s rs 


- K,(2o) [ K,(x){I _ I,(x)] a + K(x) 


+ Kolaa) [ Ksa)I2) % — tolaa) f Kz)Kaa) S- 
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The last two integrals may be evaluated as above; this time, we obtain the Wronskian 
of K, and J, and another term in 1/s*. The Wronskian term cancels the 1/zs* term in 
the above expression, leaving 


» & 


QeV> = Kolm) [| K(2)[1 — Io(a)] = + wK(e) 


cos (r — O)s — (1/s*)K,(2). 


$ sin 1s 
To evaluate the remaining integral, consider the contour integral 
© HS@[L — Jo@)\(ae/2) 


taken over the contour shown in Fig. 4. In the limit as p — 0 and R >=, the integrals 











Fic. 4. Contour for evaluation of the integral 
f 0 As(y)[1l — Lo(y)] dy/y. 


over the circular ares contribute nothing (assuming that | Re(s) | < 2; the result for s out- 
side this strip is obtained by analytic continuation). Thus, 


(2/mije***? | K(y)[1 — To(y)\(dy/y) = | H(a)[1 — Jo(x)|(dx/z). 


Equating real parts, we have 


(—2/r) sin : 78 K,(y){1 — Io(y)|(dy/y) = [ J (x)(dx/x) — | J (x) J o(x)(dx/x). 
74 “0 0 
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These integrals have been evaluated elsewhere [7, p. 391, Eq. 1; p. 403, Eq. 2], so that 


[ Kw Io(y) (dy/y) — E a 
il — Lg dy/y) = = — ft = —- oa , 
J. \¥ y Y/Y 2 sin ixs|s rs S a © 


Putting this into the expression for V, , making several cancellations, and reverting to 
the original notation, we have: 
cos (r — 6 


1y)8 ae K,(pro/c) e 


Vo = K,(pr,/0) xr ae 
apne 4s sin 37s 


2s sin 78 


a 
2. 
3 
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BOOK REVIEWS 


Approximations for digital computers. By Cecil Hastings, Jr. Princeton University 

Press, New Jersey, 1955. viii + 201 pp. $4.00. 

This is a book with a new and refreshing approach to problems of approximation. Perhaps the 
word ‘‘new”’ is misleading, since the author makes much use of the classical orthogonal polynomials of 
Chebyshev, but the newness consists in the imagination with which old methods are applied. The book 
is the outgrowth of work done by the author at the Rand Corporation in connection with its program 
of research for the United States Air Force. 

Part I, dealing with best fit, weights, iteration, solution of equations, Chebyshev polynomials, and 
other illustrative material, serves as an introduction to the set of specific approximations contained in 
the second part. The concepts are clearly and neatly presented by means of aptly chosen figures and 
graphs, while the written text is compressed to a few lines of explanation for each diagram. 

Part II consists of a list of seventy-four completely worked out approximating expressions for 
twenty-three transcendental functions of a single real variable. The approximating expressions are 
either polynimials or rational fractions except in a few cases where logarithms also are involved. Several 
different approximations are given for many of the functions, so that there are more than three times 
as Many approximating expressions as there are functions. 

The data for any given approximation are uniformly presented according to the following scheme 
(e.g. for arctan X) 


— 


Function: 





arctan X 
Range: 
-1<X <1 
Approximation: 
arctan *X = Ci\X + C3X* + C;X5 
C, = .995354 
C; = —.288679 
Cs; = .079331 


Error Curve (The graph of the error is drawn to scale, from which the error for any X can 


be seen at a glance). 
Comments: (Other approximations or pertinent information are often given under this heading.) 











In the reviewer’s opinion this book will be an invaluable piece of equipment for every large scale 


computer. 
W. E. MILNe 


Asymptotische Entwicklungen mittels der Methode der stationdren Phase. By Joachim 

Focke. Akademie-Verlag, Berlin, 1954. 48 pp. $1.14. 

The asymptotic development in 8 of the function defined by ff g(z, y)e~®*‘*'” dz dy is studied by 
the method of stationary phase. The discussion opens with a consideration of the problem of evaluating 
JS g(x)e~? ‘=) dz where, at the stationary point, g(x) behaves like (x — £)~*, 0 < \ <1. The Bessel function 
J(8) and I'(8 + 1) are discussed as examples. The double integral discussion includes consideration of 
several possibilities for the behavior of ¢ in the neighborhood of its stationary point but g(z, y) is required 
to be analytic in both z and y everywhere in the integration domain. No examples are treated for the two 


dimensional case. 
G. F. CARRIER 
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SIMPLE SOLUTIONS OF SAINT-VENANT TORSION PROBLEM 
BY USING TCHEBYCHEFF POLYNOMIALS* 
BY 


MOHAMMED M. ABBASSI 
Mathematics Department, Faculty of Engineering, Alexandria University, Egypt 


Summary. By using Tchebycheff polynomials, various forms of the stress functions 
associated with torsion problems have been obtained. The method can be applied for 
various types of cross sections. It is well known that these problems can be solved by 
other methods, but the method presented here gives a simple and direct solution in 


many cases. 
Introduction. The problem of torsion of prismatical bars by couples M, applied 
at the ends was first studied by Barré De Saint-Venant who discussed several solutions 


of the equation: 


$+ = — 2ua (1) 


ax? t ay? 
where y is a stress function, » the modulus of rigidity and a the angle of twist per unit 
length of the bar 


¥ = constant along the boundary of the cross section. (2) 


The magnitude of the torque M, on any cross section is given (See for example [1]) by 


M, = 2 |f vardy. (3) 
The only non-vanishing stresses are 22 and Yz and these are given by 
= _ - 
a-37,. #--3 (4) 


Many attempts were made to solve Eq. (1) to satisfy the boundary condition (2) for 
various cross sections, and in one of these methods the solution was obtained by the 


aid of the complex variable. 

We shall here give a simple elementary solution for Eq. (1), and the results are 
expressed in terms of Tchebycheff polynomials which are already tabulated. In putting 
y/(— Qua) = } (x* + y’) + ¢; Eq. (1) becomes 
ao , do 
aa! T oy? (5) 
Using polar coordinates, x = r cos 6, y = r sin 6, and noting that 

J 9 é sin 6 0 
Ox ~ Or r 06’ 
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*Received April 8, 1955. 
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76 
Eq. (5) reduces to the well known Laplace’s equation in polar coordinates viz: 


4 eS. Pes a (6) 


or ror ae 
An elementary solution of Eq. (6) can be obtained by putting 
@ = F,(6)-F.(r). 
It can be easily found that 
F,(0) = E, cosné + H, sin né, 
Por) = Ko’ + Lys", 


n #0 


where E, , H, , K, and L, are arbitrary constants. If n = 0, we have the solutions 
F, = dy + 1A, 
= ¢ logr+d, 


x 
| 


where dy , bo , Co and d, are arbitrary constants. The general solution of Eq. (6) can be 


put in the form 
@ = Ao + Bod + (Co + Dod) logr + >. (A,r" + B,r") cos n6 
+ (Cir" + Dr") sin né. (7) 


The stress function y is then given by 





v = z r+ A, +B,6+ (Cy + Dod) logr + > (A,r” + B,r") cos né 


—2ua t n=1 


+ (Cir + Dr") sin né. (8) 


The Tchebycheff polynomials. The Tchebycheff polynomials of the first kind 
T,(cos 6) are defined as follows: 


T,(cos 0) = cosné. (9) 
These polynomials have been treated in quite an interesting manner in [2]. We shall 
give here only what is needed in this paper. Evidently T, (cos 6) = 1 and 7, (cos 6) = 
cos 6; T,, (cos 6) for any value of n can be found by using the following recurrence relation 


T,.:(cos 8) — 2 cos 6T,,(cos 6) + T,_,(cos 6) = 0 (10) 


which can be proved as follows. 
Since cos (n + 1) 6 + cos (n — 1) 6 = 2 cos n@ cos 6 we have 


T,,.:(cos 8) — 2 cos 6T',(cos 6) + T,_,(cos 6) = 0. (11) 


By using the recurrence relation (10), we obtain the following values of T', (cos 6). 
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T,(cos 6) i. 

T (cos 6) = cos @, 

T.(cos 0) = 2 cos’ 6 — 1, 

T,(cos 0) = 4 cos’ 6 — 3 cos @, 

T (cos 0) = 8 cos* 6 — 8 cos 6+ 1, 

T;(cos 6) = 16 cos’ 6 — 20 cos’ 6+ 5 cos 8 


and so on. 
The Tchebycheff polynomials of the second kind U, (cos @) are defined as follows: 


U,(cos 6) = sin né. 


Evidently U, (cos 6) = 0, U, (cos @) = sin 6. 
Since sin (n + 1) 6 + sin (n — 1) 6 = 2 sin n@ cos 8, we obtain the recurrence relation 


U,,,,(cos 0) — 2 cos 6U,(cos 6) + U,_:(cos 6) = 0. (12) 

By using (12) we obtain the following values of U,, (cos 6) 

U,(cos 6) = 0, 

U,(cos 6) = sin @, 

U.(cos 6) = 2sin @ cos 8, 

U,(cos 6) = sin 0(4 cos’ 6 — 1), 

U,(cos 6) = sin 0(8 cos* 6 — 4 cos 8), 

U;(cos 6) = sin 0(16 cos‘ @ — 12 cos’ 6 + 1) 


and so on. 
The general solution in terms of Tchebycheff polynomials. The general solution 
(8) can be expressed in terms of Tchebycheff polynomials as follows: 


= ra 4+ Ay + Bod + (Co + Do) logr + >> (Ay + Byr")T.(cos 0) 


+ (C,r" + D,r")U,(cos 8), — (18) 


where r = (x? + y’)'” and @ = tan” y/z. 


The solutions given by By@, (Co + Do@) log r and >>*., (C,r” + D,r-") U, (cos 6) are 
in most cases not of practical importance. Putting By , Co , Do , C, and D, all equal to 
zero, we get the simple solution 


=o - if + >> (A,r" + B,r")T,(cos 8) (14) 


which is valid for cross sections symmetrical about an axis, the axis being taken as the 
initial line. 
If further we put B, = 0, we have 


= i + > A,r’T,(cos 6). (15) 


—2ya n=0 
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If the equation of the boundary can be put in either one of the forms 


i? + z. A,r"T,(cos 6) = 0 (16) 
or 
r+ DAs + Ber )T,(cos 6) = 0, (17) 


the expression on the left hand side of (16) or (17) multiplied by —2ya will give the 
stress function for the given cross section. 


1. Circular cross section. The equation of the cross section is 2? + y’ = a’ or 
1 y? — 1 q? = 0 which is of the form (16) by putting A, = —} a’ and all other constants 
equal to zero. The stress function y is given by: 
¥ = —fualr’ — a’) 


or in Cartesian coordinates 
(2 2 2 
vy = —huc(e + y — a’). 


2. Elliptic cross section. The equation of the cross section is 


ely’ 
~+4=1 
a b 


or in polar coordinates 





which can be put in the form 


Ls a’b’ ” i 
a, 3+ ae =O 


which is of the form (16) when we put A, = —4a°b’/(a’ + b’), A, = —} (a’ — b’)/ 
(a? + b’) and all other constants equal to zero. The stress function is given by: 


_ 1 2 _2a*b* : a - 0 2rn ; 
dias 1 hal rae ra Z T{c0s 0) 





or in Cartesian coordinates 


a 7; ) 
v ia Ua a” 4. b? (2 + b? . 


3. Equilateral triangular cross section. The equation of the cross section Fig. 1 
is given by 


ae - 24) 5 - 24) = 
(x + 4a)(z V3y 3 4 z+ V3 y 3 4 = Q, 


or ae ae ee , 
(x Ph ate see” aa eom ae” ind Sine = 
i.e. 5 + y) a7 4 ta (x sry) = 0 
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and in polar coordinates by 
r —-—a — —rT,(cos d=) = 0 
4 27 4a af ) 
which is of the form (16) when we put A, = —a’/27, A; = —1/4a and all” other 
constants equal to zero. 
The stress function is given by 


1 2 4 
vy=-5 tol 7 — 37 a’ — * T,(cos 0 | 


or in Cartesian coordinates 
1 3 4+, Boy ° 
y=-5 pal 2 +y - azt — a (x* — any’) | 


4. The cross section is a circle of radius a with a notch whose boundary is a circle 
of radius b whose centre is on the circumference of the circle (as shown in Fig. 2). The 
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cross section is made up of the two circles r = 2a cos 6 and r = b, its equation can be 


written in the form: 


2 2 2a 
(r° — bj) 1 — 7 C08 @ = 0 


or 
l 1 1 T.(cos lab). 
4° “a yy — 5 ar i(cos 6) + Ss T (cos 0) = 0 

which is of the form (17) when we put Ay = — 3b°, A, = —4a, B, = 3ab’ and all other 


constants equal to zero. The stress function is given by: 


1 2 2) 2a 
p= —=  — f° — —T7,(cos f 
y 5 Mar pli ; T (co | 


5. Cardioid cross section. The equation of the cross section Fig.f3 is given by 





r= (1 + cos 6) 


2a” 


Rees 


> CC 








Fig. 3. 
or 
(2a’r — cos 6)? = 1 
which can be put in the form: 
a’r’(2a’r — cos 6)? = me + cos 6) =r cos. ; 


or 


ar(2a’*r — cos 6) = r'” cos 


Nola 


which reduces to 


1 1 1 
4” - agi? Tr2(cos 6) — garTi(cos 6) = 0. 
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This equation is of the form (16) by putting A, = — 1/8a*, A, = — 1/8a’ and all other 
constants equal to zero. The stress function is given by 


l 2 I 1/2777 - 3 . 
y= “7 hal — 587 T (cos 6) — oa? rT ,(cos 0 |. 


Other forms of cross sections can be similarly obtained. 
REFERENCES 


1. S. Timoshenko, Theory of elasticity, 1934 
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BOOK REVIEWS 


Vector and tensor analysis. By Nathaniel Coburn. The MacMillan Co., New York, 1955. 
xii + 341 pp. $7.00. 


The book is divided into three parts: vector analysis (Chpts. 1-4), tensor analysis (Chpts. 5-8), 
and application of tensor analysis (Chpts. 9-12). The scope of the book is perhaps best summarized 
by giving the list of chapter headings: 1. The fundamental operations, 2. Differentiation theory, 3. 
Integration theory, 4. Vector analysis in applied mathematics, 5. Tensors in Cartesian orthogonal 
coordinates, 6. Tensors in general Cartesian coordinates, 7. Tensors in general curvilinear coordinates, 
8. The theorems of Gauss and Stokes in N-dimensional Euclidean space, 9. The differential geometry 
of surfaces, 10. Introduction to the theories of elasticity and viscous fluids. 11. The theory of compressible 
fluids, 12. The theory of homogeneous statistical turbulence. 

The reviewer feels that the book is a poor one, and this for two reasons. The first is that some of 
the author’s statements are careless and misleading. The second and more subjective reason is that the 
author’s approach to the subject is often one with which the reviewer finds himself out of sympathy. 

The first point is illustrated by the author’s statement in Sect. 8(a) that the triple vector product 
is equal to the volume of the parallelopiped spanned by the three vectors. This statement is not qualified 
by a discussion of the orientation of a vector triad or of the sign of the triple product. The author’s proof 
of the identity A. B X C = A X B. C can easily be generalized by the unwary reader to give the 
false result A. BX C=B.AX C. 

The second point is illustrated by the following examples. 

The author defines a vector as a quantity which (1) may be represented by a directed line segment 
and which (2) adds according to the parallelogram law; he carefully checks the second part of this defi- 
nition in each concrete case. This old fashioned appeal to an a priori physical concept of addition seems 
both vague and unnecessary. Thus, in the space-time of the special theory of relativity, there is no 
simple physical interpretation for the vector sum of two 4-velocities, but the vector sum of the 4-momenta 
of two particles is the total momentum of the system. In the author’s view then, the 4-momenta would be 
vectors and the 4-velocities not. Yet both 44momentum and 4-velocity transform alike under Lorentz 
transformations, and one is a scalar multiple of the other. 

The reviewer was unable to find any general discussion of contraction as a tensor operation, except 
in the special case of inner multiplication. The word “‘contraction”’ is not listed in the index. 

The author essentially restricts himself to Euclidean space of three or more dimensions and, in 
Chapter 9, to curved surfaces in ordinary space. Riemannian space is mentioned several times, but 
there seems to be no definition or explanation of the term, nor does it appear in the index. Thus the 
discussion of the Riemann curvature tensor seems rather trivial in the general context of the book, since 
the tensor reduces to a single independent component for a 2-surface and vanishes completely for a 


Euclidean N-space. 
A. ScHILp 


Machine translation of languages. Edited by William N. Locke and A. Donald Booth. 
The Technology Press of the Massachusetts Institute of Technology, Cambridge, John 
Wiley & Sons, Inc., New York, and Chapman & Hall, Ltd., London, 1955. xii + 243 


pp. $6.00. 


Although this book contains no mathematics, many mathematicians will find it to be very stimulat- 
ing and enjoyable reading. It contains a collection of fourteen essays by various authors including experts 
both in languages and in machine computers. As a result of the organization of the book there is a great 
deal of repetition of ideas but the book serves a very useful purpose in bringing these ideas to the atten- 


tion of a wider audience at a time when the subject is still in an early stage of development. 
G. F. NEWELL 


STRESS-RELAXATION IN INCOMPRESSIBLE ELASTIC MATERIALS 
AT CONSTANT DEFORMATION* 


BY 
R. 8. RIVLIN 
Brown University 


1. Introduction. The mechanics of the finite deformation of incompressible isotropic 
ideal elastic materials has received a considerable amount of attention in recent years’. 
In the present paper, we shall discuss the mechanics of the finite deformation of in- 
compressible isotropic elastic materials which behave substantially as ideal elastic 
materials when subjected to deformations during a relatively short time interval (less 
than some time 7’, say,) but, if held in that state of deformation for times much greater 
than 7’, exhibit relaxation of the deforming forces. Such a material will be called a stress- 
relaxing elastic material. It will be shown that many of the results of the theory of the 
finite deformation of incompressible isotropic ideal elastic materials apply, with slight 
modification, to stress-relaxing elastic materials. 

2. The stress in a stress-relaxing elastic material. We assume that a body of an 
isotropic stress-relaxing elastic material is subjected, at zero time, to a deformation 
in which points of the body initially at X, , in the rectangular Cartesian coordinate 
system z,; , move to x; in the same coordinate system, so that 


z; = 2,(X,), (2.1) 


where x, is independent of the time ¢, for t > 0. We assume that the stress components 
o;; in the coordinate system z, , at a point of the material initially at X, , are single- 
valued functions of the nine gradients 0x,;/dX, and of the time ¢; i.e. 


oi; = fi,(0x,/AX, , t). (2.2) 


Since the material is assumed to be isotropic in its undeformed state, the functional 
form of the dependence of o;; on dz,/0X, and ¢ must be independent of the choice of 
the rectangular Cartesian coordinate system z,; . Now, suppose x% is a rectangular 
Cartesian coordinate system related to z,; by 


zt = 4;;2; , where G; 02; = bi . (2.3) 


Then, if X* are the coordinates in the coordinate system x* of a point which is at X, 
in the coordinate system x, , we have 


X*% = a,;X; . (2.4) 
If o* are the components of the stress at time ¢ in the coordinate system xt , we have 
o* = A4,4;,0,, = f(0x%/aX*% ’ t). (2.5) 


*Manuscript received April 19, 1955. The results contained in this paper were obtained in the course 
of work on a project sponsored by the Armstrong Cork Company. The author is indebted to Dr. G. W. 
Scott, Jr. and Dr. J. T. Bergen for their encouragement and to the Armstrong Cork Company for 
permission to publish the work. 
1For reviews of the subject see C. Truesdell [1] and Green and Zerna [2]. 
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From (2.5), it can be shown, in a manner similar to that adopted in the theory of 
the finite deformation of perfectly elastic materials (see, for example, Rivlin and Ericksen 
[3, Sec. 13]), that o;; must be expressible in terms of dx,/dX, through the components 
of a symmetric Cartesian tensor, which has components C,, in the coordinate system 


x; defined by 








ee (2.6) 
ax, OX, 
and that o,;; must be a single-valued function of t and the components C,, of this tensor. 
If it is assumed that o;; is expressible as a polynomial in the nine quantities dz,/dX, , 
with coefficients which are single-valued functions of ¢, then it follows that it must be 
expressible as a polynomial in the six independent components C,, of the tensor defined 
by (2.6), with coefficients which are single-valued functions of ¢. 
If we represent the functional dependence of o;; on C,, and t by 


05; = F3(Coe ’ t), (2.7) 
then, in view of the relations (2.5), we must have 
F, (CS 9 i= 054051 F 4(C,, ; :. (2.8) 
where C;* is defined by 
* Dyk 
oo Oe (2.9) 


CS = a48alu = = 
” = Queens 
Employing the notation 6 = || ¢,; || and C = || C;; ||, Eqs. (2.7) and (2.8) imply that 
6 is an isotropic matrix function of C. 
It can be shown, in a manner similar to that adopted by Rivlin and Ericksen [8, 
Sec. 39] in the case when ¢ is an isotropic matrix function of C and the elements of 6 
do not depend on #, that a linear relation of the form 


6= aI +a,C +a,C’, (2.10) 


exists between 6, C, C* and the unit matrix I. In this relation, a, , a, and a, are single- 
valued functions of ¢ and the scalar invariants tr C, tr C’ and tr C’*. If the elements of 6 
are expressible as polynomials’ in the elements of C, then ay , a; and a, are expressible 
as polynomials in tr C, tr C’ and tr C’, with coefficients which are single-valued functions 
of ¢. 
If the material is incompressible, then the stress is undetermined to the extent of 


an arbitrary hydrostatic pressure p (say). We must then replace Eq. (2.7) by 


oi; = Fis(Coe , 2) — 184; - (2.11) 
The constancy of volume of an element in the deformation is expressible by 
det C,, = 1, (2.12) 
which leads to the relation 
det C = 3[2trC* + (tr C)* — 3trCtrC’] = 1. (2.13) 


*See also Reiner [4,5] and Prager [6]. 
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In a manner similar to that employed in arriving at Eq. (2.10) from (2.7), we obtain, 
from (2.11), 


6 = aI +a,C + aC’ — pl, (2.14) 


where a , a, and a, are single-valued functions of ¢ and the scalar invariants tr C, tr C* 
and tr C*. Since the relation (2.13) is satisfied a) , a, and a, may be expressed as single- 
valued functions of ¢ and the scalar invariants tr C and tr C*. Also, since p is arbitrary, 


we may re-write (2.14) as 
6 =a,C + aC’ — pI, (2.15) 


where a, and a, are single-valued functions of ¢ and the scalar invariants tr C and tr C’ 

If, in (2.11), F;; is a polynomial in the quantities C,, with coefficients which are 
single-valued functions of t, then, in (2.14), a> , a, and a, are polynomials in tr C and 
tr C* with coefficients which are single-valued functions of ¢. Therefore, in (2.15), a, and 
a, are polynomials in tr C and tr C’ with coefficients which are single-valued functions 
of t. 

3. Comparison with the theory of finite elastic deformations. For a perfectly 
elastic isotropic material, subjected to a deformation described by (2.1), the stress matrix 


6 is given by 


6= aI +aC +a,C’, (3.1) 


where a» , a, and a, are polynomials in the scalar invariants J, , 7, and 7; defined by 


I, = trC, I, = 4{(tr C)’ — tr C’], 
(3.2) 


I, = det C = 3[2trC* + (tr C)* — 3trCtrC’). 


a , a, and a, are therefore expressible as polynomials in tr C, tr C’ and tr C*. The co- 
efficients in these polynomials are not, in this case, functions of time as in the case of 
the stress-relaxing elastic material. They are, however, derivable from a strain-energy 
function W, which is expressible as a polynomial in J, , 7, and J; , by means of the 


relations 


_ gp aW = 2, (2% + 1,2) 
a, = 2]; al.’ <i. ai | ol, I al, 
and 
2 aw 
a. = yi al, sions 


If the perfectly elastic isotropic material is incompressible, then the stress matrix 

6 is given by 
é6=aC+a.C — pl, (3.4) 
where p is an arbitrary hydrostatic pressure and a, and a, are now expressible as poly- 


nomials in J, and J, only, while 7, = 1. Since J; = 1, it is seen from (3.2) that a, and 
a, must be expressible as polynomials in tr C and tr C’ only. a, and a, are derivable from 
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9) 


a strain-energy function W, which is expressible as a polynomial in 7, and J, only, by 
means of the relations 





ow ow ow 
= 2| — —— a = . 
ay (a } I, om) and Qe 2 al, (3.5) 


Since the stress-relaxing elastic materials considered in Sec. 2 are assumed to behave 
substantially as perfectly elastic materials, if the deformation is carried out in a suf- 
ficiently small time interval, a) , a, and a, in the expression (2.10) for the stress matrix 
must be expressible in the form (3.3) when ¢ = 0. Similarly, in the case of an incom- 
pressible stress-relaxing elastic material, a, and a, in the expression (2.14) for the stress 
matrix must be expressible in the form (3.5) when ¢ = 0. 

However, after stress relaxation has taken place, the a’s need no longer be expressible 
in terms of a strain-energy function. It is noted that, even in this case, the expressison 
for the stress matrix for the stress-relaxing and perfectly elastic materials are identical 
apart from the fact that for the stress-relaxing material the a’s depend on time and 
are not necessarily expressible in terms of a strain-energy function. Also, we shall assume, 
since this involves no added complication, that the a’s are not necessarily expressible 
as polynomials in the scalar invariants of C, but are, more generally, single-valued 
functions of these. 

4. The applied forces for a specified deformation. The body forces f; per unit 
mass of material and surface tractions F; per unit area of surface measured in the de- 
formed state, which must be applied to a body in order to maintain the stress distri- 


bution ¢;; in it, are given by 


00; ; ae 
&; + pf; = 0, 
(4.1) 
F,; - a; ;l; ] 


where p is the density of the material of the body in its deformed state, /,; are the direction- 
cosines of the normal to the surface of the body in its deformed state and it is assumed 
that the body is held in a constant state of deformation. 

Equations (4.1) apply to both the stress-relaxing and perfectly elastic materials in 
both the compressible and incompressible cases. In the theory of the finite deformation 
of perfectly elastic materials, they have been used to calculate the forces necessary 
to support a number of simple types of deformation in the incompressible case using 
the stress-deformation relation (3.4), in which a, and a, are given by (3.5). In most of 
these problems no explicit use is made of the fact that a, and a, depend on the scalar 
invariants J, and J, . In these cases we can take over the results of the theory of the 
finite deformation of perfectly elastic incompressible isotropic materials and apply 
them to the stress-relaxing elastic materials discussed in Sec. 2 by replacing dW/dI, 
and dW/dlI, by 4 (a, + I,a.) and —} a, respectively. 

We shall consider a few results in the theory of the finite deformation of perfectly 
elastic incompressible isotropic materials which, with this modification, apply also to 
stress-relaxing isotropic incompressible materials. 

(a) Simple extension. For a uniform rod of perfectly elastic incompressible iso- 
tropic material, subjected to a tensile force F, the force F is given (Rivlin [7]) in terms 
of the extension ratio \ by the relation 
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IV , 
: Law , 12m) 
sila 24( ar, * x al,/? aan 
5) where A is the cross-sectional area of the rod in its deformed state and J, and J, are 
given by 

"4 ce 1 

f- h='+y and I, = 2+ 53° (4.3) 
ix 
“a Making the substitutions } (a, + I,a2) and — } a, , where a, and a, are functions of 
S I, , J, and t, for 2W/aI, and dW/al, respectively, we obtain, for the incompressible, 


isotropic, stress-relaxing elastic material, 


P = A(n “ Lf + (x, _ |, (4.4) 
| 


where J, and 7, are still given by (4.3). Substituting for J, in (4.4), we obtain 


| F = A(n - 1) as + (x? + 1, |. (4.5) 


(b) Simultaneous simple extension and torsion of a circular rod or tube. If a tube 
of isotropic incompressible perfectly elastic material is subjected simultaneously to a 
simple extension of extension ratio \ and torsion of amount ¥, the force system which 
must be applied in order to maintain the deformation consists of (i) azimuthal surface 
tractions applied to the plane ends of the tube of magnitude © per unit area measured 
in the deformed state of the tube; (ii) normal surface tractions applied to the plane ends 
of the tube of magnitude Z per unit area measured in the deformed state of the tube; 
(iii) normal surface tractions applied to the inner curved surface of the tube of magnitude 
P per unit area measured in the deformed state of the tube. y is defined as the angle of 
twist per unit length of the rod in its extended state. 6, Z and P are given (Rivlin [8}), 
at a point distant r from the axis of the tube in its undeformed state, by 














aw <4 
= 942/24 2 i 
oa ee d al,/”: 
z — a (x2 — 1\(eW 1) | ff aw 
7* af (. \( +i a) ~ * oJt L'a (4.6) 


« «ett oe 
P= avn |r Spar. 
In these equations 7, and J, are given by 
I, =? + ; +yrr’ and I[,= u + 2. + yr’ (4.7) 


and a and b are the external and internal radii of the tube in its undeformed state. 
The azimuthal surface tractions 6 are statically equivalent to a couple M given by 


= sey [ e(H 4 180) 
M= any | (ae i al, dr (4.8) 
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and the normal surface tractions Z to a longitudinal force N given by 


pm ten — 4h) f° (2% 4 1p ce — ney [AM + 22%) a 
oka a(n 2) i (3 ‘ia °°" * i ‘Von tias@™ 9 





By taking b = 0 in the first two of Eqs. (4.6) and in Eqs. (4.7) to (4.9), we obtain the 
results appropriate to a rod of circular cross-section and radius a in its undeformed 
state. 

Replacing dW/dlI, and dW/dlI, by 4 (a, + I,a2) and — 3} a, respectively in Eqs. 
(4.6) to (4.9) and substituting the expression (4.7) for J, , we obtain the forces which 
must be applied to a tube of our stress-relaxing elastic material in order to maintain in 
it a state of simultaneous simple extension and torsion. Thus, we have 


8 


pr’ | a + (x? + : + vas | 


Z= (x? _ IE + (? + ; + v7") | + pra, 


“A/ 


+ ry [ | + (x? + : + v1) | dr 


2 


and 


aa ) 
P = —wyr | iE + (x? + x + v0). | dr, (4.10) 
while M and WN are given by 
—_ 2 l 2, 2 
M = 2ry | r E + (y a i + ~W~Ar a dr, 


N = 2n( A —_ L) [ | + (x? + 2 + v0). | dr (4.11) 


3 
— ry ae + ( + v0) | dr. 


We again obtain the results appropriate to a rod of circular cross-section by taking 
b = 0 in the expressions for 6 and Z given in (4.10) and in (4.11). 

(c) Small torsion superposed on simple extension for a uniform rod of arbitrary cross- 
section. It has been shown by Green and Shield [9]* that if a rod of isotropic incom- 
pressible perfectly elastic material, of arbitrary uniform cross-section, is subjected 
to an extension of extension ratio \ and to an infinitesimal torsion of amount y, measured 
per unit length of rod in its extended state, then the tensile force N and torsional couple 
M which must be applied in order to maintain the rod in this state of deformation are 
given by 

3See also [2]. The argument of Green and Shield given in [9] has been modified in detail by Gent 
and Rivlin [10]. 
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aw iaew 
= 9Ai— bat page 
x 24(oH +X at 4), 


(4.12) 
IW 1 aw 1 
M = 2y(2¥ 1 ew) -4a- 9} 
' Ae d dl, r (I — S)/, 
where J, and J, are given by 
2 2 1 
I,=xX + x and I,= x2 + 2. (4.13) 


A denotes the cross-sectional area of the rod in its undeformed state, J denotes the 
moment of inertia of the cross-section of the rod, in its undeformed state, about the 
axis of torsion and S denotes the geometrical torsional rigidity of the unextended rod 
when subjected to a small twist. 

From Eqs. (4.12) and (4.13), we see that for infinitesimal amounts of torsion super- 
posed on a state of simple extension, we have the following relation which is independent 
of the precise form of the strain-energy function 

r 2 
A ’ (A — 1/\))A | (4.14) 


M/y I-(U-S)/r* 

Again, following through the analysis of Green and Shield, but employing the expressions 
(2.14) for the stress in our stress-relaxing elastic material, we can obtain expressions 
identical with (4.12), in which dW/dI, and 8W/dl, are replaced by 3 (a, + J,a,) and 
— 4 a, respectively, where a, and a, are functions of J, , J, and ¢. These expressions give 
the tensile force and torsional couple which are required to support a simultaneous 
simple extension and infinitesimal torsion in our stress-relaxing elastic material. Again, 
from these expressions we can obtain the expression (4.14), which is thus seen to be 
valid also for a stress-relaxing elastic material if the extensional and torsional defor- 
mations are simultaneously applied in times small compared with those in which ap- 
preciable relaxation of stress can take place in the material, the rod is held subsequently 
in this state of deformation and the tensile force N and torsional couple M are measured 
simultaneously at an arbitrary time after the deformation is produced. 
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—NOTES— 


A NOTE ON THE SIMPLY-SUPPORTED PLATE* 


By M. STIPPES (Virginia Polytechnic Institute) 


The deflection of a simply-supported plate subjected to a concentrated load at some 
interior point is well-known in the form of a bilinear series [1]**. In a recent paper, 
B. D. Aggarwala [2] has given closed expressions for the vertical shears Q, and Q, in 
terms of Weierstrassian elliptic functions g(z) and their derivatives. The purpose of 
this note is to show how his results can be extended directly to give closed expressions 
for the bending moments M, , M, and for the resultant vertical shears V, , V, at the 
boundary of the plate. 

We begin by considering a rectangular plate the mid surface of which is oriented 
as shown in the Argand diagram of Fig. 1. If the plate is subjected to a concentrated 


e 
QO 














Fia. 1. 


load P at the point z and is simply-supported on the boundary I, then the transverse 
deflection, w(x, y), is determined by the following conditions. 


V'w = 0inR; w = 0, V’w = OonT, (1) 
and, additionally, 


a 


) — —— (2 — a)(2* — 2%) log ( — a)(e* — z%) = Biz, 20 , 2*, 2%), (2) 


167D 


where B, is analytic in R, 2* denoting the conjugate of z. 

The physical interpretation of the condition expressed by Eq. (2) is nothing more 
than the fact that the singular part of w in R is such that the vertical shear caused by 
this deflection, when summed over the lateral surface of a cylinder of radius « drawn 
about z, is equivalent to a load P. With w so determined, the bending moments, vertical 
shears, and resultant vertical shears are computed by the usual expressions 


*Received December 14, 1954; Revised manuscript received April 15, 1955. 
**Numbers in square brackets refer to the bibliography at the end of this paper. 
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M, = “Mae Hs =) M, = - (= + - =) (a) 
Q, = —D= — w Q, = -D5 Vw (b) 

(3) 
V. = —D> ie + (2 —«) a) (c) 
a - E +a-a oe. (a) 


Now in Eq. (1), we let V’w = M and find that 
V°M = 0inR; M = 0onT 


where M is subject to the condition that 


— a 
M 4nD log (z Zo) (z « o) 


is analytic in R. Now an elementary application of the Shwartz-Christoffel mapping 
theorem along with the method of images yields 





: P {lee — ele) l[ee*) — ge*))\ 
Vo = — | 2 << — te 4 
” = iD lige) — pedlle® — eedI! ad 
The function (uw) is the Weierstrassian elliptic function with periods 2w, , 2w; . (This 
result, of course, is well-known.) 
From Eq. (4), it follows that w can be written in the form 


a __ 28) log {1P@_— eedIlee*) - | 
w(x, ¥) = Ferp & — 2o)e* — 2%) I (Se, — g@)][e@ — 9%)] 
P 


T iéeD 


Bz, 20 , 2*, 2%), (5) 


where the function B is analytic inside the rectangle, R, vanishes on the boundary I, 
and satisfies the following Poisson equation in R 


a°B = | = eet) e’(2*) 








Oz dz 2+ 


oe) — 96) 9@) — 5 | + complex conjugate. (6) 

Now Eq. (4) is one linear relation between the derivatives w,, and w,, . One more 
relation between these two quantities, independent of Eq. (4), will permit their explicit 
presentation. To this end, we note that in view of the definitions of Q, and Q, along 
with Eq. (4) both Q, and Q, are harmonic functions. Consequently, there results from 
Eq. (3b) the following relation 


“0 








- -1 (2 _ &) Ff FF r 
V’°Or w) (22 oy ’ C) _ Ox" ay” (7) 
Furthermore, since 
(zt, \ _ 99 (¥ ) _ 9Q, 
V (2 0.) ~ Ox ’ V 9 Q, st oy ’ (8) 
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Eq. (7) may be rewritten as 
v{Ow + 5p (@Q. - va} = 0. (9) 


It should be remarked that this is the key relationship for the determination of the 
desired quantities. 
From Eq. (9), we have 


[yw + sp (7Q, — yQ,) = A(z, y), (10) 


where H, is a harmonic function. To determine H, , uniquely, it is necessary to know 
its boundary values, and, additionally, the nature of its singularity in R. For a con- 
venience in a future integration, we shall determine, instead of H, , the function 


} Ex : ee. =. 
T(z, y) = A,(x, y) + aDi Qg,- oD C:. (11) 
Since the plate is simply-supported at the boundary, we have 
w = 0, V*w = OonT implying [}w = OonT. (12) 
The boundary values of J, are then determined from Eq. (10) and Eq. (11). An in- 
vestigation of Eq. (5) reveals the nature of the singular part of J, . 
From the addition theorem for the g functions, it follows that g(z) = g(z*) any- 
where on I’. With this observation, the following three functions 


P 25 9"(eo) z59"(eo) + complex conjugat 
-— seceeanan 2g, Maal Os enameeime a -omplex ate 
8xD g(z*) — (2) — — 


g(z) — (Zo) 
P 9’ (Zo) 9’ (zo) . 2) 
R, = cD ee _ 2G") — pe) + complex conjugate (13) 


R, = Ry 


(considered as functions of z with z) as a parameter) vanish on I’. Additionally, each 
of these functions is a harmonic function with such type singularity in R that 


R, = 


Ji(z,y) = A(z, y) +R. + a (Q, + Rs) — SK (Q. + R) (14) 
2Di 2D 
is analytic in R. Evidently J, assumes the same boundary values on I as does J, . 
Since Green’s function for LaPlace’s equation for a rectangular region is known, 
we have 
dG 


an ds. (15) 


_ / I, 


If now we perform the indicated integration in Eq. (15), then after some laborious 


computations there results 


— : a 
F= 5 H(z, ”Y — ap (2Q, — yQ,) 
—1 a ; , 
= [4w,w3f(z) — 2(mws + wins)e — wiz*][(Q.(z, 20) — 1Q,(2, 20)] 


os 3D [4e,w3¢(2o) — 2(nws + wins)2o — wizt][Q.(2o. ,8) — 1Q,(z0 , 2)] 


+ complex conjugate, (16) 


1956] HEINZ G. HELFENSTEIN 93 


where by definition 
f@) =m, $(@s) = ms . 


(It is to be noted that ¢(z) is the Weierstrassian function associated with g(z).) 
This completes the analysis, for from Eqs. (4) and (10) it follows that 


w,, = 3V*w+ F, 
2 
Wy = 3V w — F, 
and, consequently, all of the expressions in Eqs. (3) can now be presented in closed 


form. 
If we introduce Jacobi’s theta functions, we can write 


Pi E 4(y) + rily — | 0,(v) 82(v) O3(v) Os(v) A1(v0 + v5) 8,(vo — v%) 








7 82D w, 8,(v) 6,(v + Vo) A(v = Yo) O,(v + v%) Av — vs) 
re E 41(v0) “e »| 9,(¥) 82(¥0) As(¥o) B(v0) O:(v + ¥*)0,(v — v*) 
82° D Lew, 6,(v») ee 6: (vo + v)O,(¥o — v)O,(vo + v*)O,(% — v*) 


+ complex conjugate, 


where 


Zz 20 
Pra Qu,’ ae iil 2a, 





By some arguments involving the properties of the real and complex parts of analytic 
functions z of a complex variable z, it is possible to show that 


Wey + > (yQ, + 2Q,) = ~e K,(z, 20 , 2*, 2%) + K(e , 2%), 


where K, is the harmonic conjugate of H, . However, there seems to be no immediate 
way of evaluating the function K. 
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GRAPHICAL DETERMINATION OF A DISCONTINUITY SURFACE 
BY WAVE REFLECTION* 


By HEINZ G. HELFENSTEIN (University of Alberta) 


1. Introduction. The reflection of waves, such as seismic reflection for prospecting 
has been discussed in many papers and books. The methods used in practice all share 
the simplification that the discontinuity surface is replaced in the neighborhood of the 
shot point by its tangent plane, and it is only this tangent plane which is determined. 
We shall discuss here a mathematically exact and fairly simple graphical method of 


*Received Jan. 26, 1955; revised manuscript received April 18, 1955. 
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determining points of the discontinuity surface from a limited number of well chosen 
points of the time-distance surface. The velocity of the waves in the surface layer is 
assumed to be constant and known. 

2. Notation. Introduce a set of rectangular axes in space such that the surface 
of the earth is the x, y plane and the positive z-axis points vertically downwards. Assume 
the velocity s of seismic waves to be constant and known in the surface layer down to a 
discontinuity surface D. Waves originating from a shot at the origin are reflected at D. 

To each point F = (u, v, 0) on the surface of the earth there corresponds a value T 
of the time, viz. the time from the shot to the arrival of the reflected wave at F. Multi- 
plying this time by the velocity s we obtain the total distance travelled by the wave. 
Speaking of the “‘time-distance surface” S we mean that w = sT is plotted against w, v. 

According to the law of reflection the ellipsoid of revolution which is the locus of 
the points whose sum of the distances from the two points (0, 0, 0) and (u, v, 0) is equal 
to w is tangent to the surface D. In other words: D is the envelope of all these ellipsoids; 
the surfaces S and D are related by the contact transformation associating to each point 
(u, v, w) the above mentioned ellipsoid. 

This transformation can of course be expressed by formulae connecting the equations 
of the surfaces and containing their partial derivatives. Since we can measure only a 
limited number of points of S, however, we try to avoid the computation of these partial 
derivatives. Our object is to indicate a construction by ruler and compass for points of 
D using only measurable data about S, in particular only a finite number of well chosen 
points of S. 

3. Transformation of a line-element of S. We start with the following question: 
What information about D can be obtained from the knowledge of a single (non-vertical) 
line-element of the time-distance surface, given by the two close points P = (u, v, w) 
and P, = (u + Au, v + Av, w + Aw) on a line ¢ in space? 

With each of these points an ellipsoid of revolution E and EL, is associated. The foci 
of E are 0 = (0,0, 0) and F = (u, v, 0); those of E, areO and F, = (u + Au, v + Ap, 0). 
In general, the surfaces EF and £, will intersect in a fourth order curve C, . Let J be an 
arbitrary point on C,. Then by definition of the foci, we can write 


I0+J1F =u, 
10+ I/F, =wet Aw. 


Subtracting yields 
IF, — IF = Aw. 


The last equation means that J is situated on one sheet of the hyperboloid of revolution 
H with foci F and F, and constant difference of distances from them equal to Aw. It is, 
more precisely, the sheet opening in the direction of descent of our line-element. The 
asymptotic cone belonging to this surface has its axis along the line FF, , and its apex 
is the midpoint of this segment. 

Denoting by As the distance FF, , the opening of this cone is given by cosa = Aw/As. 
On the other hand, the slope of the line PP, with respect to the u, v-plane is given by 
tan 7 = Aw/As. 

Hence, for P,; tending to P along the line ¢, the hyperboloid H degenerates in a right 
circular cone C' whose apex is at F, whose axis is the projection of the line ¢ on the u, 
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v-plane, and whose opening is found from 
cosa = tan r. (1) 


The above mentioned sheet of H containing J becomes that half of C which opens in the 
direction of descent of the line-element. 

We showed that the intersection C, of EF and £, is at the same time the intersection 
of EF and H. These two surfaces intersect in a fourth order curve K, consisting in general 
of two loops, one on either side of /, corresponding to the two halves of the cone C. The 
loop in the direction of descent of the line-element is a locus for the point R on the dis- 
continuity surface where the wave traveling to F was reflected. 

Since both of the surfaces E and C are symmetrical with respect to the u, v-plane 
the same is true for their intersection K, . Hence the projection of K, on the u, v-plane 
is a conic (more exactly, it consists of two ares of a hyperbola). 

If we know two line-elements through a point P on the time-distance surface (i.e., a 
tangent-plane) we are able to find the reflection point R by intersecting the two corre- 
sponding curves K, (both of which lie on #). All the curves K, originating from the 
o' jine-elements which make up the tangent-plane of S intersect in the same points 
(there are 2 of these points, symmetrical with respect to the u, v-plane). This follows 
from the fact that the reflection point is uniquely determined by a point of S and its 
tangent plane. 

4. Two special line-elements and the corresponding loci. For arbitrary line-elements 
this intersection of two fourth order curves is of course graphically too troublesome. 
We shall show, however, that in each tangent-plane of S there are two special line- 
elements yielding simple curves K, . 

(a) First locus for R. One of them is the line-element in radial direction: the axis 
of its corresponding cone C coincides with the axis of E, and K, therefore degenerates in 
a pair of circles the planes of which are perpendicular to the line OF. We can find them 
graphically in the following way. 

First, from the slope of the radial line-element we find the angle according to (1) by 
an obvious construction (See Fig. 1). 


ton r 











Fia. 1. 


Then assume that OF is the direction of ascent of the time-distance surface in radial 
direction. The angle @ is placed along the axis OF with its vertex at F and opening 
towards 0 (see Fig. 2). Along its other leg we mark off the distance FQ, = w. The mid- 
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perpendicular of the segment OQ, intersects FQ, in a point C, for which 
C,0 + CLF = C0; + Cr = WwW. 


Hence C, is a point of the ellipse in which the wu, v-plane intersects EZ. As the line FQ, is 
also a generating line of the cone C (because of the angle a), C, is a point of K,. Fora 
time-distance surface descending in radial direction the distance w is measured off the 
line FQ, in the opposite direction and we get other points Q, and C, . If C{ and Ci denote 
the symmetric points of C, and C, with respect to the line OF, the segments C,C! and 
CC? are the projections of the two circles of which K, consists. C,C{ is the loop con- 
taining F if S is ascending in radial direction, otherwise it is C,C3 . 

(b) Second locus for R . There is another choice of the line-element for which K, 
is still simpler. For the direction of the level-line of the time-distance surface we have 
7 = 0. From (1) we conclude therefore that a = 90°, which means that the corresponding 
cone C' degenerates in a plane perpendicular to the projection of the level-line through F. 
K, becomes then an ellipse the points of which are counted twice, and the projection of 
which on the u, v-plane is a segment which is another locus for the projection of R. This 
gives the interesting property that the projection of the reflection point is always on the 
normal n to the level-line of the time-distance surface. 

The depth of R can easily be found by means of the semi circles over C,C{ or C,C}. 

Our construction fails if C,C{ coincides with n. This can only happen if the level-line 
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is itself radial. In this case we use for instance the line-element the projection of which 
is perpendicular to OF. Its angle 7 is the same as the angle of the tangent plane of S. The 
axis of its corresponding cone C is perpendicular to OF. Hence we have to intersect this 
cone [its opening is again found from (1)] with the circle C,C{ which is done by turning 
down the plane of the latter. 

As soon as we have a point of D we have also its tangent-plane: it coincides with the 
tangent-plane of the ellipsoid LE. 

5. Choice of the points on S. The accuracy of our construction depends largely 
on the accuracy with which the directions of the level-lines and the slopes of the radial 
sections can be found. If the discontinuity surface in a certain area A is to be found it is 
therefore advisable to place the recording instruments in several straight lines L, radiating 
from the shot-point over A. From each line L, we obtain the graph of a radial section of 
the time-distance surface by joining the measured values w on L; by a smooth curve. 
From this curve the slopes r of the tangents can be found graphically. In the map of the 
area A, the level-lines can be drawn from the measured data. 

For each seismometer we obtain one reflection point with its tangent-plane. However, 
for each intersection of a level-line with one of the lines L; which occurs in a “new” 
point (where no instrument was placed) we get another point of the discontinuity surface. 
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A NOTE ON NUMERICAL DIFFERENTIATION* 
By JOHN W. MILES (University of California, Los Angeles) 


Summary. Given the matrix f = {f;}, representing f(x) at the set of points {z,;}, 
the mth derivatives of f(x) at these points are expressed in terms of all of the f; according 
to f‘” = C~'A”"Cf, where A is the sum of the skew matrix [(x,; — x;)~"] and the diagonal 
matrix formed by summing the terms in the corresponding rows of this skew matrix, 
and C is the diagonal matrix having as its elements the products of the elements in the 
corresponding rows of the skew matrix. 

1. Introduction. Let f be the column matrix 


f= {fi} = tf@)}. (1) 
We require a square matrix D such that** 
my _ a" f(a) -p" 
fe" = { a a = D"f. (2) 


*Received Feb. 28, 1955; revised manuscript received May 16, 1955. 

**The representation of derivatives in matrix form, as in (2), also has been considered by J. Kuntz- 
man in a paper presented at the International Mathematical Congress in Amsterdam (Sept. 1954), 
but no details have been published. It appears, from private correspondence with Prof. Kuntzman, 
that the results of Eq. (8) et seg. in the present paper are probably new. 
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It is common practice, especially in working with tabular data, to form derivatives 
through differencing." The usual differencing procedures, however, are not directly 
suited to the expression of a derivative at any one of a limited number of 2; in terms of 
all of the f; , as may be necessary in certain applications (e.g., the calculation of the aero- 
dynamic incidence associated with a given deflection mode in the flutter analysis of a 
highly swept or low aspect ratio wing, where accurate differentiation is important? and 
the number of control points severely limited). The desired goal of using all the f; could 
be achieved by an appropriate combination of forward and backward differencing, just as 
in interpolation,’ but it is more direct to start from the Lagrange interpolation formulaft 


fm = ef; I’ @—-2,), (3) 
C; = TI’ @; - 2)", (4) 
which expresses f(x) as the nth order polynomial prescribed by the n + 1 f, . 
2. The first derivative. Differentiating (3) with respect to x and setting x = 2; , 
we have’ 
-(#) -3 
fi ai (af —_ — p> D;j;f; ’ (5) 
D;; = Cr'C xz; — 2)", t = j, (6a) 
= >)’ (%; — 2)’; i= j, (6b) 
s=0 


where the D;; are the elements of the matrix operator D introduced in (2). Alternatively, 
we have the more symmetric form 


n 


Cit = D' (: — 2) (Cif + Cif). (7) 


i=0 


The construction of the matrix D and its higher powers is facilitated by the similarity 
transformation 


D = C’’AC, (8) 


where C, which we designate as the conditioning matrix, is diagonal and comprises the 
C,; , while A is the antisymmetric matrix formed by eliminating the factor of C;'C; from 
D,,; . If A is resolved into a skew matrix “’A and a diagonal matrix “A according to 


1—. Whittaker and A. Robinson, The calculus of observations, Blackie and Son, London, 1949, pp. 
62-68. 

tIt generally is necessary to consider at least the first and second derivatives of the wing deflection 
in calculating the wing incidence. 

*Whittaker and Robinson, loc. cit. pp. 47-50. 

ttThe prime on the product symbol, and subsequently the summation symbol, implies the omission 


of the term j = s. 
3W. E. Milne, Numerical analysis, Princeton Univ. Press, 1949, pp. 93-100. 
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A _ (A + “-. (9) 

we have the construction sequence 
A = [(1 — 8))(x, — 2)", (10) 
(A = p By 4. (11) 


C = E Il’ A. | (12) 


s=0 


We remark that, in the interests of convenience, C may be modified by the introduction 
of any constant, normalizing factor. 

3. The higher derivatives. It follows from the known properties of the similarity 
transformation, or otherwise directly from (7), that 


D” = CAC. (13) 


Few applications demand derivatives higher than the second, but we emphasize 
that, in consequence of its generation from the nth order polynomial (3), A must satisfy 


the equation 
ria _ 0. (14) 
We also remark that, in virtue of the Cayley-Hamilton theorem, (14) implies that the 


n + 1 latent roots of A all vanish. 
The elements of A* may be constructed directly from the product law 


A}? = >o A,.A,,; - (15) 


s=0 


Substituting the A,; from (10) and (11), we find 


A ie ” 24;;(A;; ‘and Ai), a ¥ d; (16a) 
=3 >’ 42, i=j, (16b) 


which affords a relatively direct construction of A’ from A. We note that the diagonal 
elements of A’ still follow from the summation of the remaining elements of a given row, 
but that A’ no longer is characterized by any direct symmetry. 

4. The case of equal intervals. We now suppose that the x; are separated by equal 


intervals of length h — viz., 
t; = to + th (17) 
so that the interpolation coefficients defined by (4) become 


C, = (-)"*/njMn — DI, (18a) 


= (—y-(" th, (18b) 
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where (;) denotes the binomial coefficient. It is expedient to remove the constant factor 
n! (— h)~" in forming the conditioning matrix, however, and we adopt the modified form 


C = | ax—(*) | (19) 


Factoring out h, the A matrix is given by 
hA = la — s)G—-—/p)' +5 do’ G-—»s } (20) 


We note that A’”’ is now characterized by the additional symmetry 


Asw.e-1 = (—)" Azz”. (21) 


In order to illustrate the foregoing results more specifically, we set n = 4 in (19) 
and (20), whence 
1 0 oO 0 0| 


25/12 —1 —1/2 1/3 —1/4] 
1 —5/6 —] —1/2 —] | 
hA = 1/2 0 _ _ >|, (23) 
1/3 1/2 5/6 —] 
1/4 1/3 1/2 25/12 | 
35/12 13/6 19/12 7/6 11/12 
-11/3 —5/3 —1/3 1/3 a 
hA)’ = 1/2 —2 —5/2 2 —1/2 (24) 
1/2 1/3 —1/3 —5/3 —11/3 | 
| 
11/12 7/6 19/12 13/6 35/12 _| 
Consider, e.g., *f(x) = sin (xx/4) and h = 1, for which f = {0, 2 + 1, 274, O}. Pre- 


multiplying this column by C~’ AC and C~'A’C, as given by (22)-(24), we find f’ = 


*This example is representative of flutter analysis, in which it often is necessary to introduce wing 
vibration modes higher (in frequency) than the fundamental and take their first and second derivatives. 

















1956] C. R. PUTNAM 101 


{0.7712, 0.5572, 0, —O. 5572, —0.7712} and f” = {0.072, —0.443, —0.614, —0.443, 
—0.072}, which differ from the exact results by {—1.80, 0.36, 0, —0.36, 1.80} % and 
{—, —1.4, 0.5, —1.4, —} % respectively. By way of comparison, we note that the 
simpler approximations to D afforded by the use of only first or first and second differences 
(forward in the first three terms and backward in the last two) yield the much poorer 
approximations f’ {0.707, 0.293, —0.293, —0.293, —0.707} and f’ = {0.914, 0.586, 
—0.086, —0.586, —0.914}, respectively. In the use of only first differences, central 
differencing in all but the first and last rows generally would be preferable—e.g., 





[1 1 0 0 0 | 
—1/2 0 —1/2 0 0 
hD = 0 —1/2 0 —1/2 0 (25) 
0 0 —1/2 0 —1/2 
0 0 0 ~f 1 | 


which yields f’ {0.707, 0.500, 0, —0.500, —0.707} for the above example. 


NOTE ON THE MANY-PARTICLE PROBLEM* 
By C. R. PUTNAM (Purdue University) 


1. It was shown by Kato [1] that the symmetric Hamiltonian operator of every 
quantum mechanical system of particles, having a potential energy of the Coulomb 
type, is essentially self-adjoint; so that its closure, H, is self-adjoint (or, hypermaximal, 
in the terminology of von Neumann [4]). In [2], Kato applied this result to prove that 
the least point, A» , of the spectrum of H for the two-electron system is an eigenvalue, 
so that A, lies in the point spectrum of H. The corresponding problem in the general 
many-electron system apparently still remains open. Thus, while the existence of a 
lower bound for the spectrum of H has been established (see [3], pp. 207-208; also [1], 
p. 205), no general criterion for the case of N = 3 electrons guaranteeing that the least 
point of the spectrum is actually an eigenvalue (or, at least, no criterion which does not 
depend upon a variational procedure with the attending convergence questions) seems 
to be known. In this connection, see [3], pp. 196-197, also the remarks of Kato [2], p. 
218, and the references cited there. 

The object of this note is to point out that \, must be an eigenvalue whenever the 
atomic number Z is sufficiently large. In fact, it will be shown that the least point, Xo , of 
the spectrum of the Hamiltonian operator, H, belonging to an atom with N electrons, atomic 
number Z, and a stationary nucleus, ts in the point spectrum whenever 


Z = 5N(N — 1)/8. (1) 


*Received March 11, 1955. This work was supported by the National Science Foundation research 
grant NSF-G481. 
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If N = 2, the inequality (1) is satisfied for all Z (where Z = 2, 3, ---); hence, as was 
already proved (and even more) in [2], \, is an eigenvalue for every two-electron system. 
If N = 3, relation (1) is seen to hold for Z = 4, thus including all three-electron ions 


except the lithium atom itself. 

A point A is said to be in the cluster spectrum of a self-adjoint operator 
H = f=. AdE(A) if either \ is an eigenvalue with an infinite multiplicity or if every 
open interval about \ contains an infinity of values (~ \) belonging to the spectrum of 
H. In order to prove the italicized assertion, it is clearly sufficient to prove that A>) < uo, 
where po denotes the least point of the cluster spectrum of H. (It is essentially this 
criterion which was verified by Kato [2] for the two-electron system.) 

2. The wave equation for the N-electron system is given by Hf = d¢ where the 
Hamiltonian 4 is defined, for a proper choice of units, by 


N 
H = Pa (-—Vi — 2Z r,) +- ba 2/T in ? (2) 
k=1 1si<kaN 


That H is self-adjoint was proved in [1]. Since the eigenvalues of each hydrogenlike 
Hamiltonian H, = — Vi — 2Z/r, are given by — Z’/n” (n = 1, 2, ---), and since the 
operator corresponding to the second summation of (2) is non-negative, it is clear that 
the least point, uo , of the cluster spectrum of H satisfies 


wo 2 —(N — 12’; (3) 


see [2]; p. 216 for the case N = 2. 

Next, let @, denote the (real-valued) normalized eigenfunction belonging to the 
ground state of the system with Hamiltonian H, and define the (normalized) function 
¢ by ¢ = [[?-: ¢ . In view of the evaluations 


2 | [(bi62)?/ra] dv; dv, = 5Z/4, (4) 
the contribution of each of the N(N — 1)/2 terms 2/r;, occurring in (2) to the inner 


product (H¢, ¢) is 5Z/4. (The integrals (4) occur in first order perturbation theory; 
see e.g., [5], pp. 163-164 for the case of the helium atom.) Consequently, 


Ao S (Hd, 6) = —NZ* + 5N(N — 1)Z/8. (5) 


If, in (5), A» = (H¢, ¢), then obviously A, is an eigenvalue. Suppose then that A, < (H¢, ¢). 
As was remarked earlier, \, is in the point spectrum of H whenever Ay < yp . In view of 
(3) and (5), this inequality is impleid by (1), so that the proof of the italicized assertion is 
now complete. 
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s A NOTE ON THE COINCIDENCE OF SOME RANDOM FUNCTIONS* 
. By GEORGE H. WEISS (Ballistic Research Laboratories, Aberdeen Proving Ground, Md.) 
Abstract. This paper deals with a system capable of randomly attaining two states 
wl “on” and “off”. For convenience this system is pictured as a pulse train. The question 
y is, what is the expected amount of time that n of such systems are in the same state 
2 simultaneously. The answer is shown to depend on the probability of a system being 
, in a given state at a specified time. A formal solution for the Laplace transform of this 
° probability is given. 
In this paper we shall consider a mathematical problem arising out of the following 
? situation: Suppose two or more new machines (relays, vacuum tubes, electric bulbs) are 


installed at time ¢ = 0, and are operated continuously thereafter. Whenever any of the 
machines break down it is replaced by a new one. Both the operating time of the machine 
and the replacement time are described by probability densities which are defined below. 
We wish to find the expected fraction of time during which all of the machines are being 
simultaneously replaced. 

A somewhat similar problem was considered’’’ in connection with the coincidence of 
periodic pulse trains coming from radar systems. 

Let us therefore consider a system which is capable of existing in two states, say 
“fon” and “‘off’’. For convenience we may picture it as a train of pulses, as shown in 
Fig. 1. As seen from the figure, the length of a single ‘‘on’’ pulse will be denoted by 
l, (the subscript k referring to the kth of n pulse trains) and the length of a single “‘off”’ 
pulse will be denoted by m, . We shall assume that the probability that an “‘on” pulse 
from the kth pulse train ends in the interval (1, , l, + dl,) is given by a,(l,)dl, . Similarly 














Fic. 1. Nomenclature for the kth pulse train. 


*Received April 5, 1955. 
1K. S. Miller and R. J. Schwartz, J. Appl. Phys. 24, 1032-1036 (1953). 
7H. D. Friedman, J. Appl. Phys. 25, 1001-1005 (1954). 
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the probability just an ‘‘off’’ pulse end in the interval (m, , m, + dm,) is B, (m,)dm, . 
Then, by definition we have 


a 2 


| ait) di « i. 
(1) 


a) 


| B.(m) dm = 1. 
70 


We may now ask, what is the expected fraction of time during which n “‘on”’ pulses, 
one from each of n trains coincide? For the answer to this question we will assume it 
possible to find the probability n,(¢) that a pulse from the kth train is on at time ¢. Then, 
in the time interval (0, 7’), the expected fraction of time during which all n pulses are 


on is 
»T n 
= | {Ty nto} dt. (2) 


Thus, the solution of the problem depends only on a knowledge of the n,(¢). The re- 
mainder of this paper gives a derivation for 7,(¢) in terms of the fundamental distributions 
a,(l,) and B,(m,). 
In order to proceed with a solution we introduce the auxiliary probability densities 
w,(t), defined by 
lim w,(¢) dt = lim Prob{pulse ends in time interval (t, ¢ + dt)}. (3) 
dt—0 dt-—0 
We shall also use the probability, A,(/), that an “‘on” pulse from the sth train is greater 


than / in duration. Hence 


A,(l) = | al) dx (4) 


I 


and A,(0) = 1, by (1 

An expression for n(¢) (we omit the subscript * in the following calculations) may be 
derived by reasoning that if a pulse is on at time ¢, then it is either the first pulse or 
a subsequent one. The probability (¢) is therefore 


t at 


n(t) = [ B(r) A(t — 7) dr + [ w(r,) dr, l B(r — 7,)A(t — 17) dr (5) 
or 
n(t) = A(t) + [ wore — r)dr (6) 
with 7 
A(t) = [ B(r)A(t — 1) dr. (7) 


70 


An expression for w(¢) is found in similar fashion. If an ‘‘on’’ pulse ends in the interval 
(t, t + dt), it must be the first of such endings or a later one, hence 


at 


w(t) = o(t) + | w(r)o(t — 17) dr, (8) 


70 
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where 


6) = | Bl)alt — 2) dr. (9) 


Since both (5) and (8) are equations of the convolution type, we may simplify results 
by using Laplace transforms. The notation will be used 








L{G()} = i] e-"'G(t) dt = G*(s). (10) 
The transform L {A(l)} is given by 
L{A(D} = += 2 (11) 
so that, from (5) and (8) 
+.) — B*@) | _1 — a*(s) | 
10 - Ol Sate ” 
and 
ree OF 1 
= a0) an - ia o ta 


which can then be substituted in (2) to give the coincidence fraction. 
An example of the use of the preceding formulae is given by a computation of the 
expected coincidence time fraction of two pulse trains, both described by 


a(l) = ae, 
(14) 


B(m) . 


where a and 6 are constants. Then a*(s) = a/(s + a) and B*(s) = b/(s + b). With these 
transforms we have 


b 


* at inna: 
so that 
@«-—— 1 - eel-e+5n (16) 
n a +. b xp / pe 


The probable fraction of time spent in coincidence is 


Cr 
- | n (t) dt 








b? { _ 2{l — exp (~(@ + BT) | 


~ (a + b» (a+ &)T 
As T —-, this approaches b’/(a + 6b)’. 
For many applications it is unnecessary to give a detailed evaluation of the integral 


[1 — exp (—2(a + 5b)T)] 
(a + BT \ (17) 
















106 NOTES [Vol. XIV, No. 1 


(2) but rather an asymptotic value as 7 —~o. In the case that each »,(¢) approaches a 
constant C, for t >, it follows that 
1 T n . 
lim =f { Il nico} dt=T[GQ. (18) 
T-@ 0 k=l k=-1 


Under certain conditions on the a,(l) and 6,(m) it is possible to show that lim ,. 
m(t) — C, , where, in fact, C, is a function of the first moments of a,(l) and 8,(m). We 
shall give sufficient conditions for 7,(¢) to approach a constant as t+ and then derive 
the value of C, . It will be shown in the following remarks that in order for n,(¢) to be 
asymptotic to a constant, it is sufficient that w,() be asymptotic to a constant. Sufficient 
conditions for w,(¢) to be asymptotic to a constant will then be quoted from the work of 


Feller’. 
At the outset let us make the following assumptions: 


(a) a(t), (tf) ~Oast-@. 
(b) The first moments, written (l,), (m,), and defined by 


(I) = [ ” Iex(I) al 
(19) 


{m) = [ mB,(m) dm 
exist. 


(c) lim A,(d = 0. 


ta 
From assumption (a) we derive 


Lemma 1: The function ¢(¢), defined in (9), approaches zero as t ~~. This lemma 
is obvious on the consideration that ¢(¢) is the probability density for the time / + m. 
With this lemma we prove 


THEOREM 1: If w(t) ~ Bast—-~-, then n(t) ~ Casto. 
Proof: Write w(t) = B + Q(t), where, by hypothesis 
| QA) | <e, e>0 (20) 
for t > N(e). The time N(e) is also chosen so that 
A(t) < «, t> N(e (21) 
which is possible because of (c). If w(é) satisfies (8), then, with a(l) and 6(m) chosen to 
satisfy (b) w(t) is a non-negative bounded function‘, hence w(t) is bounded. 
From (5) we have 
lim ¢(t) + lim “(NCE — 1)dr 


t—@ t—@ 0 


lim n(t) 


lien {B [ ade + [ Q(2)M(t — 2) ar}. 


ta 


*Feller, W., On the integral equation of renewal theory, Annals Math. Statistics 12, 243-267 (1941). 
‘This is the Theorem 2 given in the last reference. 
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The first integral in (22) can be handled by an integration by parts 


lim : (1) dr = lim {ao ~ [ ta(r) ar}, (23) 


t—@ 0 t~@ 


where the first term on the right hand side is eliminated by the use of (c). Finally 


< | [ "ann — 2) dr] + | [ “OME — 2) dr 











[ Q(r)A(t — 7) dr 











(24) 





[9 ar +e [xe- nar 


< (M + (J), 


where M is a constant independent of ¢, so that lim,.. n(t) = B¢l). 
Sufficient conditions for w(t) to be asymptotically constant have been cited by 
Feller®. Based on Tauberian theorems for Laplace transforms, these conditions are 


<e 


al) [ oat =1, 


) 


a2) / 1p(1) dt <@, (25) 


3) tim [ g(t) dt = 0. 


tw @ 


The conditions (al), (a2), and (a3) will be satisfied if assumptions (a), (b), and (c) are 
satisfied by a,(l) and 8,(m). 

The evaluation of B, makes use of the following Tauberian theorem’: If w(t) ~ B, 
as t—o then wt(s) ~ B,/s for s — 0+. The expression for «w,(s) given in (16) yields 
the result 


‘ es Soe 
wi) ~ 2 Uy) + Gm) (26) 


Therefore, if assumptions (a), (b), and (c) are satisfied on a,(l) and §,(m) an 
asymptotic result for the coincidence time fraction, (2) is 


1 af n m (,) 
al { I] mio dt~ Wa Gm 27) 


In the papers cited in Footnotes 1 and 2, the estimate (27) was not a sufficient one. 


Note added in proof: Hypothesis (c) of Eq. (19) is implied by hypothesis (b). This may 
be seen by the following steps: 


A(t) = -t | dX(zx) < [ x dv(z). (1a) 
However the integral 


- f ”  dX(2) (28) 


exists and is easily seen to equal (/), hence the last integral in (la) must approach zero 
as t approaches infinity, which proves the assertion. 


5G. Doetsch, Theorie und Anwendung der Laplace Transformation, Dover Publications, N.Y., p. 193. 
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ON DIFFUSIVE CONVECTION IN TUBES* 
By G. F. CARRIER (Harvard University) 


i. Introduction. There is a group of experiments (as exemplified by the con- 
figuration of Fig. 1) in which one determines the time dependent solute content of a 
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Fic. 1. The solute is convected into the tube at y = 0, the 
entrance condition being s = so + cos 2. 


fluid by continuously collecting the fluid through a long tube (length/radius < 1) and 
appropriately testing samples so collected for solute content. The interpretative question 
which arises is: ‘‘What is the relationship between the solute concentration at the source 
and that at the tube exit?” The answer to this question is found here for one range’ of 
values of the important parameters. The results are of interest in that the following 
unexpected’ conclusion is reached. The attenuation in the solute concentration from 
source to tube exit can be smaller for a solute with a high diffusivity than for one with a 
smaller (or zero) diffusivity. That is, the diffusivity improves the response of the “‘instru- 
ment” over the response that would be observed were the process purely convective. 
2. Analysis of the problem. The solute concentration® of the fluid in the tube is 





*Received April 23, 1955. 

‘This range is encountered in experimentation concerning the salinity of ground water in permeable 
islands. 

*T.e., unexpected by the author. 

*Actually s is the difference between the gross fluctuating concentration and the average concentra- 


tion. 
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denoted by s(r, y, ¢) and its distribution is governed by convection and diffusion ac- 
cording to the law 


8, + u(1 — r’)s, = vAs. (2.1) 


We assume here that the flow rate is not time dependent and that the Reynolds number 
is low (in the problem of footnote 1, puor./u = Re = 0 [100]) so that the flow is laminar. 
In the foregoing v is the diffusivity, wu) the maximum velocity, y and r ry are the cylindrical 
coordinates* as indicated in Fig. 1, and A is the Laplace operator. The initial condition 
states that® 

s(r, 0, 2) = exp (72). 


We wish to find the ratio of the solute discharge at y = L to the intake at y = 0. That 


is, we wish to compute 
1 1 
R= [ ra —r)se, L, 0 ar / | r(1 — r°)s(r, 0, t) dr. (2.2) 
/0 0 


It is clear that s(r, y, t) can be represented in the form 


s(r,y, ) = > F,f,(r) exp (iat — B,y) (2.3) 
and that this will be an “efficient’’ representation if Re(@,) increases so rapidly with n 
that | exp (— ByL) | « | exp (— BoL) | for moderately small N. As we shall see, there 
is a large range of the parameter w = Qrj/v for which the foregoing is true and the use 
of Eq. (2.3) is therefore to be adopted. If, then, we use Eq. (2.3) and introduce the 
notation a, = To UB,/v, Eq. (2.1) becomes 


r(rfl)’ — (iw — a,[1 — 1°] — ead) f, = 0, (2.4) 


where « = v’/uors. It can be anticipated that the results of interest (i.e., the not too 
strongly attenuated signal cases) will occur for those situations where « < 1 and it can 
also be anticipated that a will be of order w with interest centering on the case w = 0(1). 
For such eases, it is clear that the term in a; is of negligible importance and may be 
disregarded. (In the ground-water problem alluded to previously w was about 2 and e 
was 10~°.) On this basis we omit the e term with the understanding that its importance 
for any case can be checked when the answers are in by computing ea;/[iw — a,(1 — r*)). 
Consequently, we deal with the equation 


r(rf')’ — (jw — afl — r)f = 0, (2.5) 


where we have omitted subscripts. The boundary conditions require that f’(0) = f’(1) = 0 
since there must be no diffusion across the tube wall and since the origin is not a singular 
point. 

One way to obtain an accurate representation of the solutions of Eq. (2.5) for, say, 
w < 5, is to expand f in the form 


IO = Velo’, 
aa (2.6) 


with a(c) = >> ao’, 
p=0 


‘ry is the radius of the tube and r the dimensionless coordinate. 
‘We omit writing Re (real part of) for brevity. 
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where ¢ = iw. We can anticipate that the “lowest” eigenfunction will have associated 
with it an a(c) such that a(0) = 0. (Note that for o = 0; f = 1 and a = 0). Thus, the 
successive ¢, obey the equations 


-1 sy 2 9 7 
rT (re) = [p-1 — (Ayoo + Ap-191 + °** + ag,-1)(1 — r°)]. (2.7) 
Equations (2.7) can be integrated successively starting with p = O and proceeding 
monotonically to larger p. In each integration the constants of integration can be 
arbitrarily taken to be zero. (Other choices which do not violate ¢3(0) = O merely 


multiply the final result by a numerical factor.) The number a, is chosen so that the 
boundary condition ¢/(1) = 0 is satisfied. This procedure leads to the result 


fo = 1+ o(2*/8 — 27/4) + 0°(2?/96 + 52*/384 — 52°/288 + x°/256) + --- (2.8) 


a = 26 — o /24 4+ 06° /960 4+ --- 


It can be seen that the convergence is excellent for moderate values of | a | 

An alternative method of deducing Eq. (2.8) is to form the variational problem 

al 

bf {f+ [oe — all — r)If*jrdr = 0. (2.9) 
If we use the usual Rayleigh-Ritz procedure taking f to be a polynomial >> a,r” for 
0 < n < 3, and impose the proper boundary condition (f,(1) = 0), we obtain a character- 
istic equation which is cubic in a and which, in particular, gives ag = 20 — o° /24 + O(c’) 
again. It also gives a, ~ 26 + O(c) and a, ~ 160. The value for a, is certainly not to be 
trusted, but the value for a, is an excellent estimate for moderate o (say | a| < 5). If 
we adopt 26 as the order of magnitude of a, , then, we note that (6,L — 8,L), which is 
a measure of the relative importance of f, and f, , will obey the following inequality: 
Re (8, — &))L < A whenw < 24(26 — A riu)/vL). This defines, in a rough way, the 
range of utility of this analysis. The foregoing variational procedure also gives an estimate 
for the eigenfunction f,(r) which we shall not bother to record. 

The solute distribution s(r, y, t) can now be determined by noting that the ‘‘Fourier 
coefficients,” F’,, , of Eq. (2.3) are those associated with the expansion 


>, F,f.(r) = 1. 


n=0 


The orthogonality relation among the f,(r) is 


al 


| ff — r’)r dr = 0, for msn 
. 


so that F,, becomes 


al 1 
P= f.Ar\(r — r°) ar / | filr\(r — rv’) dr. (2.10) 


70 


For | ¢ | < 5 we again can show that | F, | is close to unity (we shall give an example 
soon) and’ | F, | < .01. It is clear that, for these values of o, the ratio R defined by Eq. 


*This is associated with the fact that, for | o | small, fo(r) is close to unity and all other f,(r) are 


orthogonal to fo(r). Hence f,(r) is nearly orthogonal to the function A(r) = 1. 

















111 





G. F. CARRIER 






1956] 


2.2) is dominated by errors of order 10~* exp — (8{ — 83)L (where 6’ indicated Re 8) 
by the f, contribution and is therefore given by 


= ‘al [ for(r - r°) al / ff nO — 7’) ar} exp (—a,L). (2.11) 


The evaluation of these integrals leads to the result 


— : eo 

= EE ym), a 
For the range of « we have mentioned, RF is very close to unity. In fact, for one case of 
interest, 2 = 2x/day, v = 1 em?/day, uw = 10° em/day, L = 6.10° cm, r5 = .4 cm’, so 
that ¢ = iw = 2.5i, ag = 5i — .252, Re (@L) = .038, Re (8,L) ~ 4, and it is clear that 
the second mode contribution is only 10-* exp (— 4) = 0(10~°) times that of the first 
mode. Thus 


R= 





R = .984 exp (— .038) = .946. 


This number is so close to unity that one wonders what it might have been were there 
no diffusion. The solution of Eq. (2.1) for y = O (suppressing, necessarily, the condition 
s,(1, y, t) = 0) is 


s = exp iQ{t — y/u(1 — 7°)}. (2.13) 


The efflux ratio, Ry (meaning R for v = 0) is given by 
1 
R, = 4 [ r(1 — r’) exp {—7QL/u (1 — r*)} dr |. (2.14) 


This integral may be evaluated to give 


R, (1 — mje” + m’ / u~* exp (—u) du 


= |(1 — me” — mE —m) |, (2.15) 


where m = i2L/uy . In the foregoing example, L/u. = 3/8 and R, = .91. Thus, the 
attenuation is less for the diffusive model than for that with »y = 0. However the asymp- 
totic behavior of R, with regard to m for the non-diffusing case is given by Ry ~ O(m="'). 
This implies that for sufficiently large L, the diffusive attenuation would be much 
greater than that for the non-diffusive case since the diffusive attenuation is exponential 
in, say, y. Nevertheless, it is curious that, in some low attenuation cases, the diffusive 
tr ect should decay less rapidly than the analogous non-diffusive transport. 

Aside from the foregoing result, one should also note that one can use the sampling 
technique associated with the foregoing over quite a large range of conditions and 
obtain accurate information with essentially no complicated data reduction. It appears, 
in fact, that the interpretation is so simple for these cases that a deliberate increase 
in tube length would sometimes be profitable in the experimental design. 

3. Comparison with results of G.I. Taylor. After the foregoing work was completed 
the author’s attention was drawn to the analysis of a similar problem [1]. In that paper, 
the dispersion of particular (non-oscillatory) distributions of solute was studied. In 
particular, it was noted that if a concentrated sample of solute were placed at y = 0 at 
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time zero the dispersion would occur in a manner equivalent to that in a tube with 
rigid body flow (uw independent of r) and effective diffusivity v’ = ugro/192v. If we identify 
the results of the analysis of Sec. 2 with those obtained when we write 


S. + (uo/2)S, = v’S,, , (3.1) 


(both our results and [1] imply that the transport speed is u./2), we also obtain the above 
value for »v’ provided OQr;/12v « 1. 

If now we use a Fourier synthesis over Q(i.e. if we write S = f°. H(w) >>, Fafa(w, r) 
exp (72t — 8,r) dQ) to describe the salinity distribution of Taylor’s problem, it is clear 
that the large time behavior would be governed by the small 2 behavior of the oscilla- 
tory solution and, furthermore, the only eigenfunction which would contribute appreci- 
ably would be f, . That is, the solution to the problem where a delta function distribution 
of solute was introduced at ¢ = 0 and y = 0 would be of the form 


Sy,r,) == | folr, © exp Gat — Ay) ao. (3.2) 


~ v= 


and a steepest descent evaluation of this integral would lead precisely to the result 
given by Taylor. 

The only advantages of this approach over that of [1] are that: (1) it leads to the 
observation that the equivalent diffusivity concept is valid only when the solute has 
traveled an appreciable distance down the tube; and (2) the formal procedure used here 
makes it possible to improve the accuracy of the prediction in the unlikely event that 
such improved accuracy were required. 

In so far as the problem with oscillating salinity intake is concerned, the analysis 
of [1] would have predicted the response as .962 instead of .946 since the discrepancy 


between unity and the Fourier coefficient of f, was neglected in the analysis of [1]. 
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